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AN IMPROVED MULTIPLE LINEAR REGRESSION AND DATA 
ANALYSIS COMPUTER PROGRAM PACKAGE 
by Steven M. Sidik 
Lewis Research Center 

SUMMARY 

NEWRAP is a digital computer program which can be used with ease to perform ex- 
tensive regression analyses or a simple least-squares curve fit. The program is writ- 
ten in FORTRAN IV, version 13, for the IBM 7094/7044 DCS. The major value of the 
program is the comprehensiveness of its calculations and options. 

NEWRAP computes the variance-covariance matrix of the independent variables, 
regression coefficients, t-statistics for individual tests, and analysis of variance tables 
for overall testing of regression. There is a provision for a choice of three strategies 
for the variance estimate to be used in computing t-statistics. 

Also, more than one set of responses of dependent variables can be analyzed for the 
same set of independent variables. 

A backward rejection option method based on the first dependent variable may be 
used to delete nonsignificant terms from the model. In this case, a critical significance 
level is supplied as input. The least significant independent variable is deleted and the 
regression recomputed. This process is repeated until all remaining variables have 
significantly nonzero coefficients. 

The NEWRAP program uses the triangular form of symmetric matrices throughout. 
It also allows for the use of weighted regression, computation of predicted values at any 
combination of independent variables, a table of residuals, and plots of residuals. 

By use of CRSPLT, a preregression analysis may be performed which may aid in 
the choice of model to use in NEWRAP. This program accepts the same raw data in the 
same format and computes the variance -covariance matrix and correlation matrix of all 
the variables and an eigenvector decomposition of the variance-covariance matrix cor- 
responding to the independent variables. Microfilm plots are then printed of specified 
pairs of variables. Punched output of residuals and predicted values from NEWRAP can 
also be used for more complicated residual plots than the direct use of the plotting op- 
tion NEWRAP permits. 

When a quadratic response function has been estimated (as for example in optimum- 
seeking experimentation) CREDUC may be used to obtain all information necessary for a 
canonical analysis of the function. 

The three programs together provide a useful data analysis package that can be ap- 
plied to a large variety of common research and development situations. 


INTRODUCTION 


RAPIER (ref. 1) is a very flexible multiple linear regression analysis computer 
program which has been in frequent use at the NASA Lewis Research Center. It was 
tested with the data presented in Wampler (ref. 2) and performed quite poorly. This 
alone was not very disturbing since real data are seldom even nearly as ill-conditioned 
as that set of data. A second factor, however, is that Wampler's data leads to a 5 by 5 
matrix to be inverted whereas RAPIER is designed to handle matrices of up to 60 by 60. 
With real data it is not uncommon for the matrix to become more ill-conditioned as the 
dimension increases. Often the user increases the size of the model by adding terms 
which are functions of the original independent variables (as for example in polynomial 
models) and this often leads to increased correlations and ill-conditioning. For this 
reason, RAPIER was modified primarily by rearranging the storage of variables in 
COMMON blocks and performing all the regression calculations in double precision. 

This was done without losing any of the capabilities of the original program (in fact add- 
ing new options). The resulting version is called NEWRAP. 

It may be of interest to some RAPIER users that in a number of sample calculations 
the major numerical inaccuracies arising in the regression calculations were not in- 
volved in the actual inversion of the X'X matrix but in the calculation of the inner pro- 
ducts which give 


b = (X'X) _1 (X'y) 

Thus a major improvement might be made by computing inner products in double pre- 
cision arithmetic and truncating to single precision answers without going to complete 
double precision arithmetic although the latter alternative would further increase the 
accuracy. As a matter of fact, the double precision inner product calculation is used in 
a different least-squares method proposed by Golub (ref. 3) which is reference 19 of 
Wampler's paper. 

It should be pointed out that in estimation problems an alternative to the obvious step 
of more accurate routines is provided by Hoerl and Kennard (ref. 4). They present a 
technique called "ridge regression" which uses the method of minimum mean squared 
error estimation in place of minimum variance unbiased estimation. The ridge regres- 
sion technique should have some definite appeal to statisticians, because it recognizes 
the fact that existence of ill-conditioned data indicates a problem which should be ac- 
counted for statistically as well as computationally. They do consider the problem of 
rejecting terms but their methods are not amenable to incorporation in NEWRAP in its 
present form. 

A second reason for modifying the program was the desire to provide plots of the 
residuals as was strongly recommended in chapter 3 of Draper and Smith (ref. 5). With 
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the microfilm plotting capabilities provided by CINEMATIC (ref. 6) available at the 
Lewis Research Center computer facility, this feature was also added to NEWRAP with- 
out significantly increasing printed output. CINEMATIC is a very specialized set of 
routines for the 7094/7044 DCS and 360/67 systems. If microfilm plotting is not avail- 
able at other computer installations, the subroutines used in plotting may readily be 
changed to routines which produce line printer plots or CALCOMP plots however. 

The RAPIER program used an algorithm for the coefficient calculations that inverted 
the correlation matrix and then converted this to the (X'X)~* matrix to calculate b. 

After inspection of several test cases, it seemed that this method did not improve the 
accuracy of the calculation of (X’X) - '*'. Thus it was dropped and NEWRAP inverts X'X 
directly. 

As with the RAPIER report, only the statistics and mathematics necessary to ex- 
plain the program capabilities will be presented along with illustrative input and output 
listings and listings of the programs. 


SYMBOLS 


B 

b 

b i 

b i 

b 0 

b l’ 

c 

c iJ 

D 


> b i 


E(x) 


F a, d 

V Z 1’ • • * » Z K> 
H 0 


matrix 

vector (column) 

true regression coefficient 

estimated regression coefficient 

constant term 

unknown parameters 

correlation matrix 

elements of C 

indicator variable, equal to 0 if no bg coefficient is estimated and 
equal to 1 if bg is estimated 

expected value of x (i. e. , mean of x over all possible values 
of x) 

vector of observed values minus predicted values 

statistic distributed as variance ratio with a and d degrees of 
freedom 

term of regression equation 
statistical hypothesis to be tested 
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J 

K 

k 

LOF 

M 

MS( source) 
N 

NPDEG 
N (jli, a 2 ) 

R 

REG 

REP 

RES 

r i 

S 

S c 

SSQ( source) 
S i 

TOT 

*n 

V(x) 

W, X 
w, x 

X s 

x(J) 


alternate hypothesis to be accepted if Hq is judged to be false 
number of coefficients estimated, excluding bg 
number of independent variables observed 

number of segments or cells in range of possible studentized residuals 
lack of fit 

total number of independent and dependent variables 
mean square due to source, where source is REG, RES, etc. 
number of observations 

pooled degrees of freedom for replication error 

2 

normal distribution with mean jj. and variance a 

number of sets of replicates 

regression 

replication 

residual 

number of replicates in set i 
diagonal matrix 

sum of squares correction if D = 1, and 0 if D = 0 

sum of squares due to source, where source is REG, RES, etc. 

elements of diagonal matrix 

total 

statistic distributed as Student's t with n degrees of freedom 

2 

variance of x, expected value of (x - E(x)) 

matrices 

vectors (column) 

stationary point of estimated quadratic surface 


i N 

- £ * 1 , 

N i=l 1] 


y vector (column) 

studentized residual 
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’ Z K 


e 


u 


Superscript: 


variables 

vector of observation errors 
mean of x defined as E(x) 

estimate of p. based on observation of random sample 
variance of x defined as V(x) 

2 

estimate of a based on observation of random sample 


transpose 


ESTIMATION OF BASIC LINEAR MODEL 
BASIC LINEAR MODEL 

In multiple linear regression, a dependent or response variable Y (such as temper- 
ature or pressure) measured on an object or experiment is assumed to be correlated 
with a function of one or more other variables (z^, . . . , z K ) measured on the same ob- 
ject or experiment. This function includes a number of unknown parameters 
(b^, . . bj) and can be represented as 

y = h(b l5 . . . , bj, z 1? . . . , z R ) + e (1) 

The only restriction imposed on this function is that it be linear in the parameters; that 
is, the function is of the form 


J 

y = Z Vi( z l’ • • •> Z K ) + e ( 2 ) 

j=l 3 3 1 K 

where fj(z^, • • • , Zj is a TERM of the regression equation. (A TERM is a quantity 

which may be a variable or a function of a variable, e. g. , T is a TERM and Z, after 

it is defined as Z = log T, is also a TERM. ) 

Suppose that there are N observations of the dependent variable. Let the subscript 

i indicate that the values are associated with the i n observation; in particular, the 

value of the response variable y. would depend on the observed values of the variables 

1 

( z il» - z iK^' Also, let the subscript j denote the j term in the regression 
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model so that x.. = .... z,.J) describes the transformations of the 

1] ] 11 1*- «, th 

z^, . . . , z. to produce the value of x- for the j term at the i observation. 

The regression model can now be rewritten as 


*i = Vll + Vl2 


+ Vu*‘i 


i = 1 , 


N 


( 3 ) 


where e i denotes the difference between the observed value and the expected value of 
yj. For the N observations, it is convenient to write this regression model in matrix 
notation as y = Xb + e where 



X = 



b = 





€ = 



( 4 ) 
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More often than not, the analyst feels the following model is more appropriate: 

y i = b 0 + bjX.j + • • • + bjX.j + e i i = 1, . . . , N (5) 

Let a Q = bg + bjX i + . . . + bjX j. Then, the result of adding this equation to and sub- 
tracting it from equation (5) and rearranging the terms is 

y i = (b Q + bjX J + ...+ bjX j) 

+ b^x^ - x x ) + . . . + bjtxy. -x jJ+gj i = 1, • • • , N (6) 

If, then, a dummy variable Xq is introduced such that, for all values of i, x.q = 1. 0, 
equation (6) may be written as 

y i = a o x iO + b l^ x il " x . P + • • • + b j( x u " x . j) + e i i = 1, . . . , N (6a) 

Equation (6a) now resembles equation (3) and maybe written in matrix notation, similar 
to equation (4), as y = Xb + e where now 


X = 



1.0 


1.0 



b = 



(7) 
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ESTIMATING b 


Equations (4) and (7) are similar in form and for N > J are an overdetermined set 

of linear equations. There will be some vector b which is a "best" vector to use. If 

o 

the vector e is composed of random variables such that E(ep = 0, V(e^) = cr < +°°, 
and the e i are uncorrelated, then as is well known, the method of least squares gives 
the linear minimum variance unbiased estimators b for b. And b is given by 

b = (X'X) -1 X'y (8) 

The matrix X'X divided by N - 1 is called the moment matrix of the experiment. 
The variance -covariance matrix of b is given by 

V(6) = ct^X’X)" 1 (9) 


It is important to note that when the form of equation (7) is used, X’X is given by 


X’X = 



N 


- 


E (x u - x j) 


N 

E (x u -X iKXjj -X j) 


N 


N 

E 

i 



This is seen to be symmetric and of the form 


X’X 




Hence, 
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NEWRAP uses this relation to advantage by storing only the lower triangular part 
of B and computing only the coefficients b^, . . . , bj by matrix manipulations. Then 
bQ is given by the simple equation 

b 0 = y - V 1 - b 2 x 2 - . . . - bjx j (11) 

where y = 2^/^ = a o - ^ can a * so s ^ own t* 13 -! 

V(b Q ) = V(y) + V(b'x) = + x'(X'X) 1 xja 2 

COV(b 0 ,b) = -(X'X) -1 xcj 2 


When there is no bg term in the regression model, 


N 2 N 

E x il Z x ii x i2 


X’X = 


N N 

Z x i2 x ij Z x i2 x iJ • 


N 

£ x il x iJ 


( 12 ) 


N 

z 

1 


ij 


Comparing this to equation (10) shows this form of X'X to be similar to the lower right 
submatrix in equation (10). This similarity is used to simplify notation by assuming 
that X'X represents either the form of equation (12) or the lower right portion of equa- 
tion (10) and considering the calculation of bg as a special case. Thus, further refer- 
ence to b implies 


b = 


/ b A 
Vfcj / 


There are two different methods of computing the regression coefficients which may 
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be used in NEWRAP. The first method uses bordering (ref. 7) on the full X'X matrix. 
If X'X is a nearly singular matrix, there may be problems with accuracy resulting in 
overflows or underflows causing execution to terminate without any results being printed. 
The second method uses a method of bordering which enters one term at a time into the 
model equation. After each term is entered, a full regression analysis is printed. 
Typically, if X’X is nearly singular, a number of terms will have been added to the 
model before the results become unreliable or cause execution to be terminated. Thus, 
at least a partial analysis of the full model is available to aid in selection of further 
models to submit. After all the terms have been entered, the program then switches to 
the procedure which inverts the appropriate full X'X matrix at each stage for further 
analyses. 

The use of the bordering method leads to a large volume of printed output and is not 
recommended as a standard procedure. Through use of CRSPLT as a preregression 
analysis program it may be easier to determine if bordering should be used. CRSPLT 
can also help indicate the order of arrangement of the terms of the model so that those 
thought to be most important can be entered into the model first. 

Also note that the individual observations may be weighted to perform a weighted 
regression analysis. NEWRAP permits the use of weights (ref. 5). In this case, the 
X’X and X'y matrices take the following form: 


£[ (x il - x . 1> 2 w i] 


N _ \ 

E L w i( x ii - x iX x ij- x .j)]' 


X'X = 


N p 

\ i?i L <X iJ " x . J )(x il ' x . l )w i] 


£ [Uij - X / wj 


>■ 


X'y 



e: x ii y i w i \ 



(12a) 
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CORRELATION MATRIX 


Another matrix of interest both computationally and statistically in the correlation 
matrix C. The elements of C, which are denoted Cjj, are the sample correlation coef- 
ficients between the terms and X^. These are 


N 



(x n ~ x , i>( x ij - x .j> 



(13) 


and all these numbers are between 1. 0 and -1. 0. 

The calculation of C can be expressed in matrix notation conveniently by defining 
a diagonal matrix S = diag(sj, s^, • . • , Sj) with elements 


s, =-pM= 


j = 1, 


(14) 


Then 


C = S(X’X)S (15) 

It may also be that the independent variables are random variables. Then X'X 
divided by N - 1 represents the sample variance-covariance matrix and C the sample 
correlation matrix. If the independent variables are considered to be from a multivari- 
ate distribution, it is useful in some cases to consider the eigenvalues and eigenvectors 
of X'X. For these reasons, NEWRAP includes options to compute and print these quan- 
tities. These may also be computed and printed through use of the CRSPLT program. 


ESTIMATING a 2 

For any regression model y = Xb + e, there are possibly two methods of estimat- 
o 

ing a . First, if the assumed regression model is, in reality, the true model, it is 
well known that an unbiased estimator is given by 
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_ y’y - b'x'y 
a RES(J) - N _ j_ D 

_ SSQ(RES) 
N - J - D 


= MS(RES(J)) (16) 

2 

Second, where there are replicated data points, another estimator of a , depending only 
o 

on V(ep = a for all i and not on the validity of the assumed model, is the pooled mean 

squares computed from the replicated data points. 

Assume the observations are grouped into replicate sets in sequence. Let R be 

th 

the number of sets of replicates and r^ be the number of replicates in the i Ln replicate 
set. Let 

r*+r^ 

SSQ(i) = ^ (y k - yf (17) 

k=r*+l 

where 

i-1 

r* = 53 r. 

3=1 3 

11 . __ 

It is assumed y is from the i Ln replicate set and y, is calculated only from those y 

lu n in 

in the i m replicate set. Then define the pooled svim of squares due to replication as 

R R 

SSQ(REP) = 52 SSQ(i) and the pooled degrees of freedom as NPDEG = 53 ( r i ~ !)• The 
i=l i=l 

2 

second estimator of cr becomes 


2 _ SSQ(REP) 

REP NPDEG 

= MS(REP) (18) 

It can be shown (ref. 5, p. 26) that the sums of squares due to residuals can be 
partitioned into a component due to replication and a component due to lack of fit; that is, 

SSQ(RES) = SSQ(LOF) + SSQ(REP) (19) 


12 



to use in tests of hy- 


2 

This partitioning is used later to determine the estimate of a 
potheses. 


HYPOTHESIS TESTING 
NORMALITY OF € 


As stated before, the only assumption necessary for b to be a linear minimum 

o 

variance unbiased estimator is that E(e.) = 0.0, V(e.) = a < +°°, and e { be uncorre- 
lated. If it can further be assumed that ~ N(0, cr ), a number of standard tests be- 
come available. NEWRAP computes a chi-squared statistic which can be used as an 
approximate test. 

o 

Under the hypothesis e i ~ N(0, a ), the studentized residuals defined by 





<r 


will be distributed as Student’s t with the degrees of freedom associated with the esti- 
mate a. If the degree of freedom is 30 or more, the t distribution is very close to the 
normal. 

The range of possible studentized residuals is (-°°, -h=°) and maybe divided into k 
segments or cells each with probability p i? so that each segment will have Np^ as the 
expected number of observations falling into it. Let n. denote the number of studentized 
residuals in the i tn cell. Then a chi-squared goodness-of-fit statistic may be calculated 
as 


2 

x k-l 


Z (n t - N Pi ) 2 

N Pi 

i=l 


NEWRAP computes this statistic by using an even number of cells greater than or equal 
to four and less than or equal to 20, such that the expected numbers of observations per 
cell is five or more. This statistic is not computed when there are less than 20 obser- 

i.u 

vations. The bounding values for the i in cell are Z^pZ^ where F(zp=(i-k)/N and 
F(Z) is the cumulative normal distribution function. Then each cell has the same ex- 
pected number of observations, say f = N/k. Then 


13 



2 


x k-l “ 



(Hi - f ) 2 
f 


= - E nf-N 
N i=l 1 


There is a point to be made concerning the chi-squared calculations. The validity 
of the use of the chi-squared statistic in a test depends upon the residuals forming a 
sample of independent and identically distributed random variables. This is not usually 
the case for regression residuals. Although the tail probabilities of the chi- squared 
tests might be in error, they should still be able to tell the statistician whether one in- 
tended normalizing transformation was more successful than another. 


ANALYSIS OF VARIANCE TABLE 

For most hypothesis testing of the regression model, it is convenient to summarize 
the available information in an Analysis of Variance (ANOVA) table, as follows: 


Source 

Sums of squares 

Regression 

(REG) 

Residual 

(RES) 

SSQ(REG) = b'X'y - S c a 
SSQ( RES) = y’y - b’X’y 

Total 

SSQ(TOT) = y’y - S r 


Degrees of 
freedom 


N - J - D 


N - D 


j 0 if no b n coefficient is estimated. 

S c = i — 2 u 

[Ny if a bp coefficient is estimated. 
b D = Jo if no bp coefficient is estimated. 

1 if bg is estimated. 


Mean squares 


|MS(REG) = SSQ( REG) J 
MS(RES) = SSQ( RES) (N - J - D) 


If there are replicated data points, another ANOVA table can be constructed to show the 
separation of the residual sums of squares into components from lack of fit and replica- 
tion, as in the following table: 
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Source 

Sums of squares 

Degrees of freedom 

Mean squares 

Lack of fit 
(LOF) 

SSQ( LOF) = SSQ( RES) - SSQ(REP) 

N - J - D - NPDEG 

MS(LOF) = SSQ(LOF)/(N - J - D - NPDEG) 

Replication 

(REP) 

SSQ(REP) 

NPDEG 

MS(REP) = SSQ( REP)/ NPDEG 

Residual 

(RES) 

y'y - b'X'y 

N - J - D 



CHOICE OF ESTIMATOR FOR o 2 

2 

As mentioned previously, there are two possible methods of estimating a depend- 
ing on whether there are replicated data points. This is true for any given model equa- 
tion. When the backward rejection option of NEWRAP is used, there is no longer one 

hypothetical model but a series of models. Thus, there is the choice of estimator for 
2 

cr to be made after each rejection of a term in the previous model. 

As an example, consider the model 

y = b Q + bjX^ + bgXg + bgXg + e (20) 


with replicated data points. The first step is to estimate bQ, bp bg, and bg. There 

will then be the estimators u^g^ anc * ^REP' ^ mo( * e l in equation (20) has not 

o n 

left out any important terms, UpEg^j) as well as o^gp is a valid estimator. 

The ratio F = MS(LOF)/MS(REP) can be used to test the hypothesis that there is no 
lack of fit, where F ~ F & ^ with a = N- J- D- NPDEG and d = NPDEG degrees of 
freedom. If the test accepts the hypothesis of no lack of fit, MS(RES) is a pooled esti- 
mate of a with more degrees of freedom. But there is the possibility that the hy- 
pothesis was accepted as a result of random fluctuation when there really is some lack 

~2 

of fit; that is, there is the possibility that CT jyj;g(j) is a biased estimator. If lack of fit 
is not concluded to be significant, the decision to pool or not is usually made on the 
basis of the number of degrees of freedom for replication. If this is ’’large" (no defini- 
tion of large is given herein), Ujyjp is used. If "small, ” the pooled estimate 

°RES(J) is used * 

In testing equation (20), should it be decided that bg is not significantly different 
from zero (see section t-TESTS), the coefficients of the following model would be es- 
timated: 


y = b Q + bjXj + bgXg + e 
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o 

From this model there is an estimate ^reS(J-I)' This estimate could also be biased 
since bg may be small but nonzero and the decision of bg = 0 may have been due to 
the low power of the test. 

At the first step, the lack of fit can be considered a random sample of an infinite 
possibility of biases. But the biases due to pooling mean squares after rejecting terms 
can be considered to be systematic biasing. In such a case the use of Cochran’s test 
for "the largest of a set of estimated variances as a fraction of their total” might be 
appropriate. 

NEWRAP provides three strategies of pooling estimates for use in the decision pro- 
cedure: 

(1) Never pool. This is usable only when there are replicated data points. The es- 

~ 2 

timator used in all t-tests is 

(2) Pool initial residual. This will pool the lack of fit and replication (if any) from 
the first model only. Additional mean squares due to rejected terms will be ignored. 

(3) Always pool. This strategy will always use a p E S(J-i) ^ or mo< * e l with * re- 
jected terms. 

- 2 _ 

Wherever a a or a is indicated, NEWRAP always uses the value calculated ac- 
cording to the strategy chosen by the user. 

TEST OF OVERALL REGRESSION 

One of the first tests usually applied to a regression model is the test of the overall 
significance of the model. In the notation of hypothesis testing this is stated Hq: b = 0; 
Hj: b * 0 where 



The statistic for this test is F = MS(REG)/cr . The F ~ F„ , with a = J - D, and d 

„ 2 a > a 

equals the degrees of freedom associated with a . 

Another useful statistic for judging the significance of overall regression is 

R^ = SSQ(REG)/SSQ(TOT). The sampling distribution of R does not lend itself to very 

2 

simple tests except in the case of Hq: b = 0. The main value of R is that it is a 
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number in the range 0 to 1 and 100 R is a measure of the percentage of variation in 
the y values that is accounted for by the regression model. 


t-TESTS 

In many cases, the regression model contains terms whose estimated coefficients 
are "small. " This may be an indication that the term does not have a real effect on the 
dependent variable and that the estimate is nonzero due to random sampling variation. 

If this is true, it is desirable to delete the term from the regression model. A test 
statistic for deciding this is 


t = 



( 21 ) 


— t +U __ 1 

where (X'X)^ denotes the i Ln diagonal element of the (X'X) - matrix. The statistic 
t ~ An equivalent test statistic is 


2 b i 

F = r = (22) 

a 2 (X*X)- 1 

where F ~ This * s °^en referred to (ref. 5) as the partial F-test. The 

quantity b i/[(X ,x )i/] is called the additional sum of squares due to bp if were last 
to enter the equation. NEWRAP computes and prints the t-statistics, the probability as- 
sociated with the interval (-t, t), and the additional sums of squares for each term. 

This particular test is the basis for the rejection option of NEWRAP. The analyst 

^ Q 

initially chooses which a estimator to use by the choice of strategy. Then the analyst 
may choose a confidence level which all coefficients must meet. For example, suppose 
a confidence level of 0. 900 is chosen. The t-statistic is then computed for each coeffi- 
cient, and the coefficient with minimum |t| is identified. If min|t| > q 950’ 

terms are concluded to be significant at the 0. 1 level (or 90. 0 percent level of confi- 
dence). If min|t| < b N-J-D 0 950’ the term corresponding to the minimum |t| is 
dropped from the hypothetical model, and the regression is recomputed. This process 
is repeated until all remaining coefficients are significant at the specified level of prob- 
ability. This procedure can be overriden by an option which allows certain specified 
terms of the model to be retained regardless of the significance of the coefficient. Ken- 
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nedy and Bancroft (ref. 8) present a study indicating the backward deletion method is 
slightly more efficient than forward selection in special situations. 

PREDICTING VALUES FROM ESTIMATED REGRESSION EQUATION 

Regression equations are often used to predict an estimated response at some con- 
dition of the independent variables. Useful estimates of parameters to know are the 
variance of the regression equation and the variance of a single further observation at 
the desired combination of the independent variables. 

Let x’ = (x,, . . . , x T ) denote the vector of independent variables at which a pre- 

± 7 J o 

diction is desired. Let x* = x - x. Let o _ denote the estimated variance of the 

~2 M- x 

regression equation at x. Let <7 x denote the estimated variance of a single further 
observation at x. Then, 


' } M-x = 52 g + X *'( X ’ X )" lx *] 

Jl.O +- + x+^X’X ) -1 x*] 


%2 -2 
a - a 
Y-x 


(23) 


(24) 


where, as before, D = 1 if a bp coefficient is estimated and D = 0 if a bQ coefficient 
is not estimated. The quantity s = is called the standard error of estimate and 

often is used as a simple approximation to o x . This approximation is close if N is 
very large and x = x, in which case, 





s 


2 


When x * x, this may be a poor approximation. NEWRAP accepts input vectors x and 

- ~ -2 ~ ~ 2 

computes y = b Q +b 1 x 1 + . . . + bjXj, as well as x , cr^ x , cr y x , a y> x , and the 

standard error of estimate. 
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N EWRAP PROGRAM 
USERS GUIDE TO NEWRAP INPUT 
Illustrative Regression Problem Requiring No Transformations 

The illustrative example is described in chapter 7 of reference 5. The data is re- 
produced in table I. Figure 1 presents this data in a sample input form. 

The basic model to be fitted is 


y = b Q + b lXl + b 2 x 2 + b 3 x 3 + b 4 x 4 


(25) 


TABLE I. - OBSERVED VALUES 


FOR EXAMPLE PROBLEM 


Unit 

number 

X 1 

x 2 

*3 

x 4 

y 

1 

-75 

0 

0 

-6 b 

1. 4 

2 

175 

0 

0 

l;vi 

26. 3 

7 

0 

0 

-65 


29. 4 

8 

0 

0 

165 

-6 

9.7 

9 

0 

0 

0 

150 

32. 9 

10 

-75 

-75 

0 

150 

26. 4 

11 

175 

175 

0 

-65 

8. 4 

14 

-7 5 

-75 

-65 

150 

28. 4 

15 

175 

175 

165 

-65 

11. 5 

18 

0 

0 

-65 

-65 

1. 3 

19 

0 

0 

165 

150 

21. 4 

20 

0 

-75 

-65 

-65 

. 4 

21 

0 

175 

165 

150 

22. 9 

24 

0 

0 

0 

-65 

3.7 

3 

0 

-75 

0 

150 

26. 5 

5 

0 

-75 

0 

150 

23. 4 

16 

0 

-75 

0 

150 

26. 5 

4 

0 

175 

0 

-65 

5. 8 

6 

0 

175 

0 

-65 

7. 4 

17 

0 

175 

0 

-65 

5. 8 

12 

0 

-75 

-65 

150 

28.8 

22 

0 

-75 

-65 

150 

26. 4 

13 

0 

175 

165 

-65 

11. 8 

23 

0 

175 

165 

-65 

11. 4 
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SAMPLE INPUT 


FORTRAN STATEMENT 


13 SAMPLE NEWRAP PROBLEM 2 


DATA IS FROM DRAPER AND SMITH. | [ | . 

APPLIED REGRESSION [ANALYSIS (REFERENCE 4 OF NEWRAP REPORT) 

C-HA PTER. 7 I . . ! .1 

INITIAL MODEL EIQUATIOM IS 1 

Y = jBO + (B1X1 + ,B2X2 .+, B3X3. + B 4 X.4 . .+ jERR. j 

I i • ■ ■ ■ * 

Y = CHAMBER PRESSURE I j 

XI = TEMPERATURE OF CYCLE. I 

X 2. o ,V.I BRA T.I ON, LEVEL 

, , . X 3. =; .DROP ( SHOCK,), . j , . . j 

X4 = STATIC FIRE j j 

RESIDUALSl ARE. BEING REQUESTED .TO Jb.E. PUNCHED. ,( F.O.R, .CRjSPL.T. PROGRAM) 

FOR RESIDUAL PLOTTING ANALYSES I ' 


L . . . 1 . 1 
, 5F6 . 0 ) 


- 65 . , 1.. 4 

, .1.5,0 , ,2,6,. 3 

. .1,5.0 , ,2.9... 4 

-65 _9 . . . 

5«40 41 4?|45 44 4? 46 4J 46 49 50 Si Si 


SAMPLE INPUT (Concluded) 


FORTRAN STATEMENT 


„ » 

»»*»» vm 9J 

SI iZ Si 54 JS J* 

0 

0 

0 

150 

7,5 

,-7.5 

,0 

, .1.5.0 

75 

, 17.5 

c 

-.6.5 

J,5 

-7.5 

- 65 

, ,15.0 

75 

17.5 

1 65 

,-.65 

0 

0 

-.6 5 

-6.5 

0 

0 

1 6.5 

1.50 

0 

-7.5 

- 65 

-65 

0 

175 

165 

150 

0 

0 

0 

-65 

0 

-75 

0 

. i 1 , 5 , 0 

,0 

-7 5 

0 

, , i.s.c 

0 

- 75 

0 

150 

0 

175 

o 

- 65 

,0 

, i 1 , 7 , 5 

p 

. -.6.5 

_0 

.. . .1.7? . . 

, ,c 

, . .-.6,5 

0 

. . .-.7.5 , . 

.-6.5 

, . ,1.5,0 

,0 

- 75 

,-,6,5 

15.0 

_0 

1 75 

165 

, ,-65 

_0 

175 

16 5 

. ,“, 6 . 5 


Figure L - Sample input form 


T 



The preceding model requires no transformations of the tabulated data for the dependent 
or independent variables. Suitable input statements are also given in figure 1. 

A subsequent example will illustrate the requirements on the input cards when trans- 
formations are involved. 


Detailed Description of Input Cards 

This section of the report describes the input cards as classified into nine sets ac- 
cording to table II. 


TABLE n. - FUNCTIONS OF INPUT SETS 


Set 

number 

Name of set 

Purpose 

1 

IDENTIFICATION 

Identify and describe problem 

2 

PROBLEM SIZE 

Define problem size 

3 

LOGIC 

Specify general logical controls 

4 

MODEL 

(a) MODEL SIZE 

(b) TERMS 

(c) TRANSFORMATIONS 

(d) CONSTANTS 

Define model equation 

5 

REJECTION 

Backward rejection controls 

6 

REPLICATES 

Identify replicated data 

7 

FORMAT 

Give data format 

8 

DATA 

Input observed data 

9 

PREDICTIONS 

Predicted values data 


The model equation is defined by set 4 of table II. An example illustrating the use 
of one blank card for input set 4 which can be used for simple linear regressions is pre- 
sented by figure 1 and table I. A second example illustrating the use of the set of 
MODEL cards in the presence of prior constants and transformations will be given at 
the end of this section of the report. A pictorial representation of an input deck is given 
by figure 2. 
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DATA FOR PREDICTIONS 
(FMT) 



DATA 

(Input set 8) 


FORMAT 
(Input set 7) 



FMT 

(13A6) 


REPLICATION 
(Input set 6) 



NARAY 

(1914) 


REJECTION 
(Input set 5) 


^DELETE 

P 

(LI) 

(F3.3) 


Figure 2. - Concluded. 



Details of the input cards are as follows: 

(1) IDE NTIFIC ATION (I, IDENT)(I2, 13A6): IDENT is Hollerith data used to identify 

the problem. I indicates the number of additional cards to read for further 
identification or description (columns 1 to 78). 

(2) PROBLEM SIZE (NOVAR, NODEP, NOTERM, NOOB, NTKEEP)(3I4, 15, 14) 


NOVAR 

NODEP 

NOTERM 

NOOB 

NTKEEP 


Number of input independent variables (number of z's in eq. (2)) 
Number of input dependent variables 

Number of terms in model equation (number of x's in eq. (3)). 
Note that bg is not counted as a term. 

Number of observations 

First NTKEEP independent terms of model equation will be re- 
tained in model regardless of significance level 


(3) LOGIC : One card with nine one-column fields 

BZERO b Q term appears in model equation (T or F) 

IFTT t-statistics and their descriptive confidence levels are to be com- 

puted (T or F) 

IFWT Weight of 1. 0 is applied to all observations (T or F). If this is 

F, see input sets 7 and 8 for further information. 

IFCHI Compute and plot residuals (T or F) 

STORYX Calculate eigenvalues and eigenvectors of X'X (T or F) 

IFSSR Model shall be increased by one term at a time using bordering 

method for matrix inversion (T or F) 

ECONMY Use economy version of output (T or F). NEWRAP does not 
print X’y, (X'X)”*, or C when set to F. 


ISTRAT Pooling strategy is 1, 2, or 3: 

1. Never pool. Use replication error as estimate of error. 

If 1 is selected and no replication is found, strategy 3 
is used. 

2. Pool initial residual only. 

3. Pool all residuals. 
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PNCH Punch residuals and predicted values (T or F). If T, observation 

number is punched and then residuals and predicted values are 
punched in (16, 4E16. 8/(6X, 4E16. 8)) format in pairs (observa- 
tion number, e^, y^, eg, Yg, etc. ). 

(4) MODEL : The MODEL cards are used to manipulate the observed input data, 

supplied by input set 8, into the form of the desired model equation. There are 
four subsets of this input set 4, namely, MODEL SIZE, TERMS, TRANSFOR- 
MATIONS, and CONSTANTS, of which the latter three are used only in the de- 
velopment of complex models. 

If a simple linear model is being analyzed, the MODEL SIZE card is left blank, in- 
dicating that the number of transformations is zero and the number of constants to be 
read in is zero. In this case, the TERMS, TRANSFORMATIONS, and CONSTANTS 
cards of this input set are not expected by the program, and the program assumes the 
independent and dependent variables are arranged on the input cards of input set 8 as 

x l» x 2’ * ' • » X J’ y l» • * *’ y NODEP 

where NODEP is the number of dependent variables. 

If a weighting factor other than 1. 0 is to be used (i. e. , if item 3 of the LOGIC card 
contains an F), the value of the weighting factor for each observation must appear as the 
last item in the list, so that in this case the data for each observation is entered on the 
cards as 


x l» x 2> ■ • •> x j> y !>•••> y NODEP’ WT 

If the weighting factor is identically 1. 0, NEWRAP reads a total of M numerical values 
for each observation, where M is the sum of the number of independent and dependent 
variables. The variables are stored consecutively in an array called X, beginning with 
location 01 and ending with location M. If the weighting factor is not identically 1. 0, 
then M + 1 numerical values are read for each observation, but the last value, being 
the weighting factor, is treated and stored separately. The data in X are used with 
their appropriate weighting factors to cumulatively create X'X and X'y as shown in 
equations (12a). 

The remaining discussion of this set explains the use of transformations and/or con- 
stants to build more complex models. Therefore, the reader who does not immediately 
need a complex model may skip this material and proceed directly to the description of 
input set 5. 

As mentioned previously, there are up to four subsets of the MODEL cards. Their 
purpose is to give the structure of the model equation and thereby specify the initial 
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operations to be performed on the input data. As used here, CONSTANTS means any 
numerical value specified to be in the model equation in advance of parameter estimation. 
These numerical values are read from the CONSTANTS cards. 

Also, the word TRANSFORMATIONS is to be interpreted as the operations per- 
formed on the input data (read from data cards) to compute the fj values (eq. (2)) of the 
model equation. The structure of these functions (and of any transformations of the de- 
pendent variables) is read from the TRANSFORMATION cards. Finally, the word 
TERMS is to be interpreted as the computed results of the operations specified by the 
transformation (including any operations that leave the input data unchanged). The re- 
sults of the TRANSFORMATIONS are stored in an array CON, and the TERMS cards 
designate the order of the relative locations in CON where the final values for the terms 
of the model equations are to be found. 

The four subsets, MODEL SIZE, TERMS, TRANSFORMATIONS, and CONSTANTS, 
will be described in detail now. Also, at the end of the description of this input set, a 
summary of these cards, with the formats used, is given for convenience. 

The MODEL SIZE card specifies NTRANS and KONNO(2I4) where 

NTRANS Number of transformations that will be performed 

KONNO Number of constants that will be read in which are required to specify model 
equation 

If the number of transformations is zero, and therefore, the number of constants is 
zero, the TERMS, TRANSFORMATIONS, and CONSTANTS cards are not expected by 
the program. This being a simple linear model case, only the MODEL SIZE card, 
which can be blank, is necessary in this subset, but the values for the observations 
which are provided in input set 8 must conform to the order as specified in the first 
three paragraphs describing this input set. 

When, however, a more complex model is desired, information must be supplied 
instructing the program as to (1) where to find the values for the TERMS of the equation, 
(2) how to create the terms from the variables and the constants, and (3) what the values 
of the constants are. This information is supplied on the TERMS, TRANSFORMATIONS, 
and CONSTANTS cards. 

The numerical values to be used in the transformations are stored in two arrays 
called X and CON. The transformations always require that an operator (some value 
from CON) performs an operation (see table m) on an operand (some specified value 
from X) to produce a result which will be stored in CON. Thus CON serves two pur- 
poses. First, if the number of constants (KONNO) specified on the MODEL SIZE card 
is nonzero, that many constants will be read from the CONSTANTS cards and stored in 
CON beginning with location 01. If the number of constants is zero, a CONSTANTS 
card is not expected by the program. Second, all intermediate and final results of 
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TABLE III. - OPERATIONS* AND CODE NUMBERS 


[X indicates a value from X and C a 
value from CON.] 


Operation 

Resulting 

Operation 

Resulting 

code 

operation 

code 

operation 

(OP) 


(OP) 


00 

No operation 

16 

1. 0., SQRT(X) 

01 

N + C 

17 

C**X 

02 

x*c 

18 

10. 0**X 

03 

c/x 

19 

SINH(X) 

04 

EXP(X) 

20 

COSH(X) 

05 

X**C 

21 

(l.O-COS(X)), 2.0 

06 

ALOG(X) 

22 

ATAN(X) 

07 

ALOGIO(X) 

23 

ATAN2(X, C) 

08 

SIN(X) 

24 

X**2 

09 

COS(X) 

25 

X**3 

10 

SIN(«*C*X) 

26 

ARC SIN( SQRT(X) ) 

11 

COS(77*C*X) 

27 

2. 0*!7*X 

12 

1. 0/X 

28 

1. 0, (2. 0* ;r*X) 

13 

EXP(C/X) 

29 

ERF(X) 

14 

EXP(C/X** 2) 

30 

GAMMA(X) 

15 

SQRT(X) 

31 

X C 


a AU function names and operations are consistent with 
FORTRAN IV mathematical subroutines. 


transformations are also stored in CON as specified on the TRANSFORMATION cards. 
The TERMS card then specifies which of the locations in CON finally contain the values 
needed to construct the X'X and X'y matrices. After all the transformations have been 
performed on an observation, the contents of the relative locations of the CON array 
specified on the TERMS card are moved back to X in consecutive locations beginning with 
location 01. 

Note especially that CONSTANTS data are stored in CON from location 01 through 
KONNO. Thus, if a transformation specifies that a result is to be placed in any of these 
locations, the result will replace the constant, so that further operations on subsequent 
transformations would use the new value stored instead of the constant value to which it 
was initialized. Care should be taken, therefore, that the results of the transformations 
be stored in relative locations greater than KONNO. 

Each transformation code is made up of four subfields of two card columns each, 
with the following interpretation: 
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Subfield 

Interpretation 

1 

Operand 

Relative location in X 

2 

Operation (OP) 

Arithmetic operation 

3 

Operator 

Relative location in CON 

4 

Result 

Relative location in CON 


Thus, subfield 1 always references the X array, and subfields 3 and 4 reference the 
CON array. The result of every transformation is a term which is stored in the desig- 
nated location of the CON array, with the added feature that, if the term is stored in 
relative location 61 or beyond, it is also stored in the parallel location in the X array. 
This is illustrated by the arrows in figure 3. This feature allows successive transfor- 
mations to be performed more easily. 

The OP (operation codes) are tabulated in table HI. The transformation with 
OP = 00 is simply an identity transformation. This transfers data from X to CON so 
that when terms are selected there is a value available in CON that can be moved back 
to X. 

When there are no transformations, NEWRAP assumes the first NOTERM values 
on a data card are the independent variables and the last NODEP values are the depend- 


Location X 


Location _ CON 


01 

02 

03 

04 

05 

06 


56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 
67 


95 

96 

97 

98 

99 



Input constants 
followed Dy terms 
for regression 
equation 


Data stored are 
stored in same 
relative location 
in both arrays 
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Figure 3. - Map of X and CON arrays. Data transferred into any location of CON 
array beyond location 60 are immediately duplicated in same relative location in 
X array. 





ent variables. When transformations are used, this convention need not hold for the raw 
input data but instead holds for the terms on the TERMS card. Thus, the first NOTERM 
values input on the TERMS card indicate the locations of the CON array which corre- 
spond to the independent variables and the last NODEP values indicate the locations of 
the CON array which correspond to the dependent variables. If the analyst desires to 
force certain terms to remain in the model regardless of their significance, these terms 
must be the first terms of the model. Then if the input for NTKEEP of the PROBLEM 
SIZE card is not zero, the first NTKEEP terms of the model will be retained. 

The complete sequence of MODEL cards and the formats used are summarized as 
follows: 

(a) MODEL SIZE (NTRANS, KONNO)(2I4): NTRANS specifies the number of trans- 
formations required and KONNO the number of CONSTANTS involved. NTRANS may not 
be greater than 100 and KONNO not greater than 60. If NTRANS = 0, the following three 
sets are skipped and the program goes directly to input set 5. 

(b) TERMS (4012): One or more cards as necessary, using two-column fields to 
denote the relative locations of the CON array containing the final values for the terms 
to be used and the order in which they enter into the model equation. The number of 
terms used is specified on the PROBLEM SIZE card. 

(c) T RAN SFORMATION S (4012): As many cards as necessary containing up to 

10 transformation instructions per card. Each transformation instruction is composed 
of four two-column fields. See table in for the list of available transformations. 

(d) CONSTANTS (5E15. 7): As many cards as necessary containing the number of 
CONSTANTS as specified by KONNO. Up to 60 CONSTANTS may be specified. If 
KONNO = 0, these cards are not expected by the program. 

(5) REJECTION (DELETE, P)(Ll, F3. 3): If DELETE is set to T, the backward re- 
jection option is used and the desired level of confidence is given by P. The P value is 
written without a decimal point so that a 95 percent confidence level is indicated by a 
950, a 99. 9 percent level as 999, and so forth. 

(6) REPLICATION (REPS)(L1): If REPS is F, the program skips to set 7 and as- 
sumes there is no replicated data. If REPS is set to T, then more cards are read in 
2014 format specifying: 

IREP in the first field of the first card indicates the number of replicate sets. 

NARAY in the second field of the first card and the remaining fields of this and 

succeeding cards consists of an array containing the number of observa- 
tions in each of the replicate sets. 

Note that it is not safe for the program to assume that all the data points for an ex- 
periment with the same levels of the independent variables are true replicates. Thus 
the user must explicitly specify the truly replicated sets. NEWRAP does check that all 
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independent terms within a replicate set are the same. If not, the program stops. A 
nonreplicated data point is considered to be a group of size 1. Note that the data in ta- 
ble I are grouped to clearly indicate the replicated data points. 

(7) FORMAT (INPUT, FMT)(I2, 13A6): INPUT specifies the unit number on which the 
input data is stored; and FMT supplies the format for reading it. 

Note that, if a weighting factor other than 1. 0 is to be used, its value will be read 
with each data point, and the format must allow for this. 

The example from Draper and Smith (ref. 5) uses a weighting value of 1. 0 for all 
data. The format is (12X, 5F6. 0) since there are four independent and one dependent 
variable to be read. If a weighting value other than 1. 0 is used, it must appear with 
every data point as the last value on the card. In such a case, the format could, for 
example, be (12X, 5F6. 0, F10. 3). 

(8) DATA : Each observation consisting of the given z's and y's read by the exe- 
cution of one READ statement. Thus, there will be at least one card for each observa- 
tion. As mentioned previously, if the transformation option is not used, the program 
expects the first variables to be the independent variables, in the order in which they en- 
ter the model, followed by the dependent variables and then the weighting value if 
IFWT = . F. Otherwise, if transformations are used, the independent and dependent 
variables may be entered in any convenient order, because the TERMS card(s) will be 
needed to specify the order in which the values will enter the model equation. However, 
if IFWT = . F. , the weighting value is still the last value supplied with each observation. 

(9) PREDICTIONS (PREDCT)(L1): 

(DATA) If, with the program LOGIC (input set 3) card, a 

computation of residuals is requested by a T in card column 4, then predicted values of 
the dependent variables are computed for all the input values of the independent vari- 
ables. In addition to these predicted values, predictions at other values of the indepen- 
dent variables might be desired. In PREDICTIONS input, one card with one column is 
used to indicate if these other predictions are desired (T or F). If this is F, a new case 
is started and the new case should start with input set 1 cards. If it is T, the following 
cards are read: One card with one four-column field specifying the number of predic- 
tions desired. This is followed by cards with the values of the independent variables 
at which predictions are desired. Only the final regression model is used, but the num- 
ber of independent and dependent variables originally supplied on the PROBLEM SIZE 
data cards are read. All transformations indicated on the MODEL cards are performed. 
Then the proper terms are chosen by the program to correspond to the final model. 

Since the dependent variables are not needed in this part of the program, the numerical 
values for the dependent variables are dummy values and should be in the appropriate 
range so that when subroutines required for the transformations (e. g. , ALOG, SQRT) 
use these values, abnormal exits will not occur. 
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illustrative Problem Requiring Transformations 


As an example of the MODEL cards usage consider the following. Suppose the 
model we are required to construct is 

+ 273 - 15 ) = b o + b l z l + b 2 z 2 + b 3 z 3 + b 4 z l z 2 + b 5 z l z 3 + b 6 z 2 z 3 

+ b 7 Zj + bgZg + b 9 z| 


Thus, in terms of equation (2) we have 


X 1 

x 2 

x 3 

x 4 

x 5 

x 6 

x ? 

x 8 

x 9 

y 


= z, 


= z, 


= Zr 


- ZjZ 2 


- Z lZ g 


- z 2 z 3 


= z 


= Z f 


= Z r 


"3 

= iog^y + 273 - 1& ) 


Table IV shows a sequence of transformations which could be used to construct this 
model equation. Figure 4 shows how the MODEL cards describing this equation would 
appear on a FORTRAN data sheet. Figure 5 shows the X and CON array contents both 
before and after the transformations are performed upon one observation and the X ar- 
ray after the appropriate terms have been selected according to the TERMS card data. 
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TABLE IV. - SEQUENCY OF TRANSFORMATIONS FOR EXAMPLE 



Interpretation 

x 1 - CON(ll) 

x x - X(61), CON(61) 

x 2 - CON( 12) 

x 2 - X(62), CON(62) 

x 3 - CON( 13) 

Xg - X(63), CON(63) 

x^ - CON( 17) 
x 2 - CON( 18) 

x^ - CON(19) 

x lX g - CON( 14) 

x l x 3 "* CON(15) 

x 2 Xg - CON( 16) 

y + 273. 15 - X(98), CON(98) 

log 1Q (y + 273. 15) - CON(20) 




Figure 4. - An example of MODEL cards. 
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SAMPLE NFWRAP PROBLEM 

DATA IS FROM DRAPER AND SMITH 

APPLIED REGRESSION ANALYSIS (REFERENCE 4 OF NEWRAP REPORT) 
CHAPTER 7 

INITIAL MODEL EQUATION IS 

y = chamber pressure 

Xl= TEMPERATURE OF CYCLE 
X;= VIBRATION LEVEL 
Xj = DROPISHOCK) 

X4= STATIC FIRE 

Y - BO + B1X1 + B2X2 + B3X3 + B4X4 ♦ ERR 
RESIDUALS ARE BEING REQUESTED TO BE PUNCHED (FOR CRSPLT PROGRAM) 

FOR RESI3UAL PLOTTING ANALYSES 


C» •*•••••••*•••• *•*•••****•*• 

4 1 V 24 

TTTTTFF1T 

THERE IS A BO TO ESTIMATE 
THERE ARE 18 REPLICATE SETS 

111111111111113322 


(12X.5F6.0) 

SAMPLE NCWRAP PROBLEM 

TERMS OF THE EQUATION, OBSERVATION = 1 

-75.00000 0 0 -65.00000 1,400000 

*» REPLICATE SET 1 *********************•*****************•*< 


TERMS OF THE EQUATION, OBSERVATION = 2 

175.0000 G 0 X50.0000 26.30000 

** REPLICATE SET 2 .*...*..*..*••******•*•**•***************< 


TERMS OF THE EQUATION, OBSERVATION - 3 

0 0 -65.00000 150.0000 29.40000 

** REPLICATE SET 3 ***************************************** 


TERMS OF THE EQUATION, OBSERVATION = 4 

0 0 165.0000 -65.00000 9.700000 

** REPLICATE SET 4 ***************************************** 


TERMS OF THE EQUATION, OBSERVATION = 5 

000 150.0000 32.90000 


** REPLICATE SET 5 


TERMS OF THE EQUATION, OBSERVATION = 6 

-75.00000 -75.00000 0 150.0000 26.40000 

** REPLICATE SET 6 

TERMS OF THE EQUATION, OBSERVATION = 7 

175. 0003 175.0000 0 -65.00000 8.400000 

•* REPLICATE SET 7 ***** .................................... 

1ERMS OF THE EQUATION, OBSERVATION = 8 

-75.00000 -75.00000 -65.00000 150.0000 28.40000 

** REPLICATE SET B ***************•*****. ****•**•••**•***•**< 


TERMS OF THE EQUATION, OBSERVATION = 9 

175.0000 175.0000 165.0000 -65.00000 11.50000 

»* REPLICATE SET 9 ************* »•••••••«»•••••*•••••••• »**«< 

TERMS OF THE EQUATION, OBSERVATION = 10 

0 C -65.00000 -65.00000 1.300000 

** REPLICATE SET 10 ..*.***. •••••••••»••••••••»•»•••••••**•••< 

TERMS OF THE EQUATION, OBSERVATION = 11 

0 0 165.0000 150.0000 21.40000 

** REPLICATE SET 11 **..****.*****.***.....*.****•********»*.« 


TERMS OF THE EQUATION, OBSERVATION = 12 

0 -75.00000 -65.00000 -65.00000 0.400000 

*• REPLICATE SET 1 c *•»*••******•******»*#*.********»•**»*»•*( 


TERMS OF THE EQUATION. OBSERVATION = 13 

0 175.0000 165.0000 150.0000 22.90000 

** REPLICATE SET 13 » •••••••••••••••••»•• •••••*•••••••••• ••••< 


TERMS OF THE EQUATIUN, OBSERVATION = 14 

0 C 


-65.00000 


3.700000 



*• REPLICATE SET 


l* 


ItRMS OF THE FQUATIUN, OBSERVATION = 15 

0 -75.00000 0 150.0000 26.50000 

TERMS OF THE EQUATION, OBSERVATION * 16 

0 -75.00000 150.0000 23.40000 

TERMS OF THE EQUATION, OBSERVATION = 17 

0 -75.00000 C 150.0000 26.50000 

*• REPLICATE SET 15 

DEP. VAR. 1 SSQ= 6.4066836 SUM* 76.400000 MEAN= 25.466666 


TERMS 

0 

OF 

THE 

EQUATION, OBSERVATION * 
1 75.0000 0 

18 

-65.00000 

5. 800000 

ItRMS 

0 

OF 

THE 

EQUATION, OBSERVATION * 
175.0000 0 

19 

-65.00000 

7.400000 

TtRMS 

0 

OF 

THE 

EQUATION, OBSERVATION = 
* 75.0000 

20 

-63.00000 

5. 800000 


• * REPLICATE SET 16 

DEP. VAR. 2 SSQ= 1.7066677 SUM= 19.000000 ME AN= 5.3333333 


TERMS OF THE EQUATION, OBSERVATION « 21 

0 -75.00000 -65.00000 150.0000 2B. 80000 

TERMS OF THE EQUATION, OBSERVATION = 22 

0 -75.00000 -65.00000 150.0000 26.40000 

** REPLICATE SET 17 

DEP. VAR. 1 SSQ* 2.8800044 SUM= 55.200000 MEAN* 27.600000 


TERMS OF THE EQUATION, OBSERVATION * 23 

0 1/5.0000 165.0000 -65.00000 11.80000 

TERMS OF THE EQUATION, OBSERVATION = 2 4 

0 175.0000 L65.0000 -65.00000 11.40000 

•• REPLICATE SET 18 

DEP. VAR. 1 SSQ= 0. 8000 137E-01 SUM= 23.200000 MEAN* 11.600000 


MEANS OF IN3EP AND OEP VARIABLES 

12. 5000 ^0 33.333333 25.000000 42.500000 16.579167 


X TRANSPOSE X MATRIX WHERE X IS DEVIATION FROM MEAN 


ROW 

1 

105000.0 




ROW 

2 

62500.00 

263333.3 



ROW 

3 

26250.00 

115000.0 

173700.0 


ROW 

4 

-26875.00 

-161250.0 

-49450. 00 

277350.0 


X TRANSPOSE V MATRIX WHERE X AND Y ARE DEVIATIONS FROM MEAN 


RUW 

1 

-1103.750 

ROW 

2 

-12398.33 

ROW 

3 

-2767.500 

ROW 

4 

25875.25 


CORRELAT ION 

COEFFIC IfcNTS 




KOW 1 

1.000000 




ROW 2 

0.375865 

1.000000 



ROW 3 

0.194372 

0.537706 

l.COOOOO 


ROW 4 

-0.157485 

-0.596668 

-0. 225295 

1.000000 


THE FOLLOWING ARE EIGENVALUES OF X TRANSPOSE X MATRIX 

489039.32 169G38.49 99358.063 81947.348 


THI S 

IS THE 

MODAL MATRIX DR 

MATRIX OF EIGENVECTORS. 

EIGENVECTORS ARE WRITTEN IN COLUMNS LEFT TO RIGHT IN SAM CROER AS 

1 IGENVALUES 

ROW 

1 

0. 177713 

0.201943 

0.820685 

0. 504097 

ROW 

2 

0.674835 

0 .220435 

0.212133 

-0.671670 


ROW 

3 

0.368164 

0.666753 

-0. 500091 

0.420793 


ROW 

4 

-0. 6*0269 

0.682691 

0. 177159 

-0.343240 
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SAMPLE <=‘WRAP PKUBLfM 

EACH COLUMN CONTAINS THE COEFFICIENTS FOR ONE DEPENDENT TERM 
CONSTANT TERM ( BO I 
II. o *123 

REGRESSION COEFFICIENTS (Bl,... f 8K) 

1 0.7iU17E-C2 

2 0 • 1 1 41356-0 1 

3 0.42ri5S3t-02 

4 G.lUi.374 


ANOVA OF REGRESSION ON DEPENDENT VARIABLE I 


SOU (CS 

SUMS OF SOUARFS 

DEGREES OF FREEDOM 

MEAN SQUARES 

REGRr >S I ON 

2462. 56296 

4 

613. 64 739 

RESIDUA. 

249.0 366- 8 

19 

13. 107.898 

TOTAL 

2711.59958 

23 



R SQUARED = SSOtREG) / SSQITOTI = 0.908139 R = .952974 

STANDARD ERRUR OF tSTIMAT. 3.620385 

USING POOLING STRATEGY I THF .-RRJR Mt AN SQUARE * 1.0455595 WITH DEGREES OF FREFOOM = 6 

F = MSI :EG)/MS(tRR)= 333.5* COMPARE TO F{ 4, o) 


ANOVA OF LACK OF FIT 


SOURCE 

SUMS OF SQUARES 

DEGREES OF FREEDOM 

MEAN SQUARES 

LACK 3F FIT 

237.963249 

13 

18.3048649 

REPLICATION 

11.0733571 

6 

1.84555951 

RESIDUAL 

249.0366 08 

19 

13. 1071B98 

F = Ms(LOF)/MS(REPS) = 9.918 


SUMS OF SQUARES DUE TO EACH VARIABLE IF IT WERE LAST TO ENTER REGRESSION 

1 5.20 '730 

2 14.2*010 

3 2.22 ’ 554 

4 1783.654 


STANDARD DEVIATION OF REGRESSION COEFFICIENTS (DERIVED FROM DIAGONAL ELEMENTS OF (X TRANSPOSE XIINVERSE FAT K IX ) 

0 0.354538 

1 O.454375E-02 

2 0.406457E-02 

3 0.392621E-02 

4 O.32c>OS0E-O2 


(X TRANSPOSE X) INVERSE MATRIX 
ROW I 0.111867 E— 04 

ROW 2 -0. 320414E-05 0.895205E-05 

ROW 3 0.22021 5E-06 -0 . 426 580t - 05 C. 826766E-05 

ROW 4 -0. 739629E-06 0.413364E-05 -0. 904698 E-06 0.576159E-05 

SAMPLE N c WR AP PROBLEM 

CALCULATE; T STATISTICS 

THE T STATISTICS CAN bE USED TP TEST THE NET REGRESSION COEFFICIENTS Bill. 

1.679485 
2.782620 
1 . 09 71 4o 
31.0878V 

UNDER NULL RYPOTHESIS THE INTERVAL (-T,TJ WHERE T IS GIVEN ABOVE, HAS APPROX PROBABILITY LI5TFC BELOw. 
MINUS SIGN INDICATES PROB EXCEEDS .999. 

1 0.856 

2 0.968 

3 0.685 

4 -0.999 


THE DESIKfD VALUE OF PROBABILITY IS 95.0 PERCENT 
HE TERM X ( 3) IS BLlNG DELETED 
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SAMPLE .J-rf*AP PSUHLt-M 

EACH COLUMN CONTAINS THE COEFFICIENTS FOR ONE DEPENDENT TERM 
CONSTANT TERM (BO) 

11.7 439 

REGRESSION COEFFICIENTS (Bl,...,CiK) 

1 0.751732E-02 

2 G.13 217E-01 

4 0.10.834 


ANOVA OF REPRESS IUN ON DEPENDENT VARIABLE 1 


SOURCE 

SUMS OF SQUARES 

DEGREFS OF FREEDOM 

MEAN SQUARES 

REGRESS! ON 

246D. 34143 

3 

820.113808 

RESIDUAL 

251.256162 

20 

12. 5629079 

TOTAL 

2711.59958 

23 



R SQUARED = SSQ( REG) / SSQI TOT) = 0.907339 R = .952544 

STANDARD ERROR OF ESTIMATE 3.544419 

USING POOL INS STRATEGY 1 THE ERROR MEAN SQUARE = 1.8455595 WITH DEGREES OF FREEDOM = 6 

F=MS(KES)/MS< ERR)= 444.37 COMPARE TO FI 3» 6) 


ANOVA OF LACK OF FIT 


SOURCE 

SUMS OF SQUARES 

DEGREES OF FREEDOM 

MEAN SQUARES 

LACK OF FIT 

240. 184803 

14 

17. 1560571 

REPLICATION 

11.0733571 

6 

1.84555951 

RE S I DUAL 

251.258162 

20 

12.5629079 

F = MS(L3F)/MS(RfcPS) = 9.296 


SUMS OF SQJARES DUE TO EACH VARIABLE IF IT WERE LAST TO ENTER REGRESSION 

1 5.053833 

2 27.08253 

4 1839.095 


STANDARD DEVIATION OF REGRESSION COEFFICIENTS (DERIVED FROM DIAGONAL ELEMENTS OF (X TRANSPOSE X UNVERSE MATRIX ) 

0 0. 35 1 33 4 

1 0.45*25 5E-02 

2 G. 35 £.93 3E-02 

4 0.327752E-02 


(X TRANSP.JSE X) INVERSE MATRIX 

ROW i 0.1 1 1808F-04 

ROW 2 -0.3U9U52E-05 0.675106E-05 

ROW 3 -0,713401 E— 06 0.362557E-05 0.564431^-05 

SAMPLE N-WRAP PROBLEM 

CALCULATEu T STATISTICS 

THE T STATISTICS CAN BE USED TO TEST THE NET REGRESSION COEFFICIENTS oil). 

1.654799 

3.930729 

31.5673b 

UNDER NULL HYPOTHESIS THE INTERVAL (-T»T) WHERt T IS GIVEN ABOVE, HA 5 APPROX PROBABILITY LISTEC BELOh. 
MINUS SIGN INDICATES PROB E XC F EDS .^99. 

1 3.851 

2 3.991 

4 -3.999 

THE DESIRiD VALUE OF PROBABILITY IS *5.C PERCENT 

THE TERM a( 1) IS BEING OELtTtD 
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SAMPLE 4 • WRAP PROBLEM 

EACH COLU'-’N CONTAINS THE COEFFICIENTS FOR ONE OE PE NO ENT TFRM 
CONSTANT rERM IB') 

11,7 871 

REGRESSION COEFFICIENTS (Bl,...,BKl 
2 0.155995E-01 

4 0.102354 


ANOVA OF REGRESSION ON DEPENDENT VARIABLE 


SOU <CE 

SUMS OF SQUARES 

DEGREES OF FREEDOM 

MEAN SQUARES 

REGRESSION 

2455.28763 

2 

1227.64381 

RESIDUAL 

256. 311962 

21 

12.2053314 

TOTAL 

2711.59958 

23 



R SUUARED = SSQ(REG) / SSQITOT! = 0.905476 R = .951565 

STANDARD ERROR OF ESTIMATE 3.493613 

USING POOLING STRATEGY 1 THE ERROR MEAN SQUARE = 1.8455595 WITH DEGREES OF FREEDOM = 6 

F=MS( REG)/MS(£RR)= 665.19 COMPARE TO FI 2 , 6) 


ANOVA OF LACK OF FIT 


SOURCE 

SUMS OF SQUARES 

DEGREES OF FREEDOM 

MEAN SQUARES 

LACK UF FIT 

245.238604 

15 

16.3492401 

REPUCAT ION 

11.0733571 

6 

1.84555951 

RESIDUAL 

256. 311962 

21 

12.2053314 

F = MS(L0F)/MS(REPS1 = 8.859 


SUMS OF SQUARES DUE TO EACH VARIABLE IF IT WERE LAST TO ENTER REGRESSION 
2 41.26723 

4 1871.546 


STANDARD DEVIATION OF REGRESSION COEFFICIENTS IDERIVED FROM DIAGONAL ELEMENTS OF (X TRANSPOSE X J INVERSE MATRIX J 
0 0.354375 

2 0.329B93E-02 

4 0. 321 44BE-02 


(X TRANSPOSE X) INVERSE MATRIX 
ROW 1 0 • 589681 E-05 

ROW 2 0.34283b£-05 0.559879E-05 

SAMPLE NEWRAP PROBLEM 

CALCULATED T STATISTICS 

THE T STATISTICS CAN BE USED TO TEST THE NET REGRESSION COEFFICIENTS BID. 

4. 728664 
31.84463 

UNDER NULL HYPOTHESIS THE INTERVAL t — T * T I WHERE T IS GIVEN ABOVE » HAS APPROX PROBABILITY LISTED BELOW. 
MINUS SIGN INDICATES PROB EXCEEDS .999. 

2 D.997 

4 -0.999 

THE DESIRED VALUE OF PROBABILITY IS 95.0 PERCENT 
SAMPLE NEWRAP PROBLEM 

FOR EACH DEPENDENT TERM AND OBSERVATION IS PRINTED 
OBSERVED RESPONSE (Y OBSERVED! 

CALCULATED RESPONSE tY CALC! 

RESIDUAL ( YDBS- YCALC=YD IF 1 
STUOENTIZLD RESIDUAL ( Z ) 

Y OBSERVED 1.4000 

Y CALC 5.C350 

Y DIF -3.6550 

STUOENT I ZED -2.6905 

Y OBSERVED 26.300 

Y CALC 27.C63 

Y DIF -0. 7o33 

STUOENTIZlD -0.5O1V 


Y OBSERVED 29.400 

Y CALC 27.063 

Y OIF 2.3367 

STUDENT I ZrD 1.7200 
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Y OBSERVED 9.7^00 

Y CALC 5.';S5U 


Y DIF 

STUDENT l £-0 

4.&450 

3.4191 

Y OBSERVED 

Y CALC 

Y DIF 

STUDENT I ZED 

32.900 
27.^63 
5. 8367 
4.2964 

Y OBSERVED 

Y CALC 

Y OIF 

STU06NTI Z*-D 

26.400 

25.893 

0.5J66 

0.3729 

Y OBSERVED 

Y CALC 

Y DIF 

STUDENT I Z C D 

B. 4000 
7.78^0 
0.6*50 
0.4527 

Y OBSERVED 

Y CALC 

Y DIF 

STUDENTI Z‘ » 

26.400 
25.093 
2.5*66 
1 .8*51 

Y OBSERVED 

Y CALC 

Y DIF 

STUDENT 12* 0 

11.500 

7.7B50 

3.7150 

2.7346 

Y OBSERVED 

Y CALC 

Y DIF 

S TUDENT I ZCD 

l.ioOO 

5.C550 

-3.7550 

-2.7641 


SAMPLE NE HR 

AP PROBLEM 

Y OBSERVED 

o 

o 

Y CALC 

27.063 

Y DIF 

-5.6633 

STUDENTI ZED 

-4. 1688 

Y OBSERVED 

0 . 4000 

Y CALC 

3 • 8u51 

Y DIF 

-3.4851 

STUDENT IZt-D 

-2.5654 

Y OBSERVED 

22.900 

Y CALC 

29.793 

Y DIF 

-6.8932 

STUDENTI ZED 

-5.0741 

Y OBSERVED 

3. 7000 

Y CALC 

5.0550 

Y OIF 

-1.3550 

STUDENT I Z ED 

-0.9974 

Y OBSERVED 

26.500 

Y CALC 

25.893 

Y DIF 

0.6066 

STUDENT IZED 

0.4466 

Y OBSERVED 

23.400 

Y CALC 

25.893 

Y DIF 

-2.4934 

STUDE NT I Z L D 

-1.8354 

Y OBSERVED 

26.500 

Y CALC 

25.893 

Y DIF 

0.6066 

STUDENTI Z.-D 

0.4466 

Y OBSERVED 

5.8000 

Y CALC 

7.7850 

Y OIF 

-1 . 9a 50 

STUDENTI /CD 

-1.4611 

Y OBSERVE) 

7.4C00 

Y CALC 

7.7 8 5u 

Y OIF 

-0. id 50 

STUDENTIZ; D 

-0.2834 

Y OBSERVfcn 

5 .8000 

Y CALC 

7.7850 

Y DIF 

-1.9850 

STUDENTI ZlO 

-1.4611 


SAMPLE NlHRAP PROBLEM 

Y OBSERVED 28. BCG 

Y CALC 25.893 

Y I) t F 2.9o66 

S TUOtNT 1 Z. n 2.1396 



Y OBSERVED 

26.400 

Y CALC 

25.893 

Y OIF 

0.5066 

STUDENTIZED 

0.3/29 

Y OBSERVED 

11.800 

Y CALC 

7.7850 

Y OIF 

4.0150 

STUDENT I ZtO 

2.9555 

Y OBSERVED 

11.400 

Y CALC 

7.7850 

Y OIF 

3.6150 

STUDENTI ZED 

2.6610 


CHI -SQUARE STATISTIC WITH 1 DEGREES OF FREEOOM = 5.000000 

SKEWNESS = 8.087309 

KURTOSIS = 76.73972 


HISTOGRAM I 


FREQUENCY 

9 

2 

5 

8 


* 





• 



• 

' HrJH 

• 



• 


• 



• 


• 


• 

* 


• 


• 

• 


• 



* 


• 

* 

• 

• 


• 

• 

• 

» 

INTERVAL 

1 

2 

3 

4 


CLASS 

*01* UNIT05, EOF. REC= 00000 FIL = 00002 


NEWRAP DOCUMENTATION AND LISTINGS 


The contents of this section include a flow chart of the program, a listing of the 
routines used in NEWRAP and their major functions, the call structure of the program, 
a dictionary of the program, and the listing. 


General Mathematical and Logical Flow of Program 


The flow of operation in NEWRAP is illustrated in figure 6. 
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Routines and Their Major Functions 


FORTRAN 

name 


Function of routine 


BORD 

EIGEN 

HIST 

INVXTX 

LOC 


MATINV 

MFIX 

NEWRAP 

OUTPLT 

PREDCT 

RECT 

RSTATS 

SUMUPS 

TRAN 

TRIANG 

TTEST 


Inverts symmetric matrix of order n by adding bordering column to 
already inverted matrix of order n - 1 

Computes eigenvalues and eigenvectors of input symmetric matrix 

Prints histogram of residuals 

Inverts symmetric matrix 

When given row and column indices of symmetric matrix element, it com- 
putes location this element would have if only lower triangular part were 
stored as vector. 

Controls inversion process; computes regression coefficients; computes 
eigenvalues and eigenvectors of X*X if requested 

Prints X'X and computes and prints C 

Executes overall problem control; computes replication error; controls 
deletion of variables when given results of t-test 

Computes residuals at observed points and plots them. Compute chi- 
squared statistic 

Computes predicted values, variances, and standard deviations of regres- 
sion line and further observations at specified points 

Writes rectangular matrix 

Computes regression statistics and writes regression and lack-of-fit analy- 
sis of variance tables 

Constructs X’X and X'y matrices one observation at a time, in double 
precision 

Performs transformations 

Writes lower triangular part of symmetric matrix 

Computes t- statistics and their significance levels; determines which vari- 
able should be deleted 


Call Structure of Program 

The call structure of the program is illustrated in figure 7. 
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j pRIANG 

Figure 7. -Call structure of NEWRAP. 


Dictionary of Program 


FORTRAN 

name 


Mathematical 

symbol 

Description 

ACTDEV 

e i 


Error in observation i; difference between observed 
and predicted response 

B 

b 


Regression coefficients other than the constant 

BO 

b 0 


Constant regression coefficient 

BZERO 



Logical variable set to T if constant bg coefficient 
should be in regression model 

CHISQ 

2 

X 


Chi- squared statistic 

CON 



Constants used in transformations, and results of 
transformations 

DELETE 



Logical variable set to T when deletion of terms is 
desired 

DUMMY 



Extra array used for plotting data 

ECONMY 



Logical variable indicating suppress printout of 
X’X, X'X deviations, and C if T 

ERRMS 

S 2 


o 

Estimate of a used in hypothesis tests 

FMT 



Variable input format 
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FMTTRI 

IDENT 

EDOUT 

IFCHI 

IFSSR 

IFTT 

IFWT 

INPUT 

INPUT5 

INTER 

IOUT 

JCOL 

KONNO 

LENGTH 

LIST 

NARAY r-j 

NCON 

NERROR 
NLOF n 

nodep 

NOOB N 


Format for printing matrix 

First identification printed at top of each page 

Original sequence number of each term relating re- 
duced models to original model 

Logical variable set to T if residual computations and 
plots are desired 

Logical variable set to T if sequential regressions 
are desired 

Logical variable set to T if t-statistics are desired 

Logical variable set to T if all weights of observa- 
tions are 1. 0 

Input logical tape unit number for data 

Set equal to 5 to denote input device is card reader 

Tape unit where input data is stored for use in 
OUTPLT 

Sequence number of term among those remaining 
which is to be deleted 

Total number of independent and dependent terms in 
regression model 

Number of constants originally supplied for transfor- 
mations 

Number of locations in matrix storage area currently 
needed 

Set equal to 6 to denote output device is printer 

Number of replications per replicate set 

Array containing addresses in CON array for use in 
transformations 

Degrees of freedom for error mean square estimate 

- J - NPDEG - D Degrees of freedom for estimating variance due to 

lack of fit 

Number of dependent variables 

Number of observations 
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NOTERM 

J 

Number of terms in current regression model 

NOVAR 

K 

Number of independent variables to be read 

NPDEG 

NPDEG 

Pooled degrees of freedom for replication error 

NREG 

J 

Degrees of freedom for determining variance due to 
regression 

NRES 

N - J - D 

Degrees of freedom for estimation of residual vari- 
ance 

NTERM 


Array containing locations of terms in CON array 
that should be in regression model 

NTOT 

N - D 

Total degrees of freedom 

NT RAN 


Array containing transformation codes for use in 
performing transformations 

NTRANS 


Number of transformations to perform 

NWHERE 


Location in X array of first dependent variable; used 
in prediction routine to adjust for deleted terms 

NXCOD 


Array containing addresses of variables (or terms 
with address >60) for use in transformations 

P 


Probability that interval (-t, t) must have before a 
term is considered to be significant 

PNCH 


Logical variable set to T if residuals are to be 
punched 

POOLED 

SSQ(REP) 

Array containing pooled sums of squares from rep- 
lications for each dependent term 

PREDCT 


Logical variable set to T if prediction option is de- 
sired 

REPS 


Logical variable set to T if there are replicate sets 
in data 

REP VAR 


Array containing replication variance of each de- 
pendent term 

RESMS 


Array containing residual mean square or variance 
of each dependent term 

RNLOF 


Reciprocal of degrees of freedom for lack of fit 

RNREG 


Reciprocal of degrees of freedom for regression 
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RNRES 


Reciprocal of degrees of freedom for residual 

RWT 


Reciprocal of total weight 

SATRTD 


Logical variable indicating that there are no degrees 
of freedom for residual if T 

STORYX 


Logical variable set to T if eigenvectors and eigen- 
values of X’X are to be computed and printed 

SUMX 

2>,2> 

Array containing sums of independent and dependent 
terms 

SUMX2 

E* 2 ,£y 2 

Array containing sum of squared independent and de- 
pendent terms 

SUMXX 

X’X 

Sums of squares and crossproducts matrix, and 
variance-covariance matrix of independent terms 

SUMXXI 

(X’X) -1 

Inverse of variance -covariance matrix of independent 
variables 

SUMXY 

X’y 

Array containing sums of crossproducts of independ- 
ent terms with dependent terms 

TOTWT 

w i 

Sum of weight of observations 

X 


Before transformations are performed, this contains 
the variables as read in. After transformations 
are performed, appropriate data from CON array 
are placed here according to information on 
TERMS cards. 

XCHK 


Array used in checking if all values of independent 
terms are the same within a replicate set 

ZEAN 

E(X), ii 

Expected or mean value (or X) 


Program Listing 

SIBFTC BLDV 

BLOCK DATA 1 

COMMOM /FRMTS/ FMT < 1 3 > » FMTTR I ( ) 

COMM0')/BIG/B(60»9)»SUMXY ( 60 » 9 ) .SUMXX (1830) » SUMXX I ( 10 30) 

X, DUMMY (2300) 

DOUBLE PRECISION B , SUMXY .SUMXX , SUMXX I 

COMMO'J/MED/BOIO) .SUMX (69 ) , SUMX 2 ( 69 ) , SUM Y2 < 9 ) ,ZtAN(69) , 

X CON( 99 ) ,ER RMS (9), I DENT (13) , I DOUT < bO > , NCON ( 200 ) » NTERM ( 69 ) , 
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X NTRANI 100) .NXCODI 100 ) , POOLED! 9) .REPVARI 9) ,RESMS( 9) ,X( 99) b 

DOUBLE PRECISION 80,SUMX,S0 MX2 , SUMY2 , ZE AN 9 


COMMON/SMALL/ 

BYPASS, BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

10 

X IFTT, 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

11 

X ISTRAT, 

JCOL, 

KONNO, 

LENGTH, 

LIST, 

12 

X NERROR, 

NODEP, 

NOOB , 


NOTERM, 

13 

X NOVAR, 

NP DEG, 

NRES , 

NTRANS, 

NWHERE , 

1A 

X P, 

PREDCT , 

REPS, 

RWT, 


15 

X STORY I , 

STORYC, 

STORYX, 

TOTWT, 

WEIGHT, 

16 

X ERRFXD, ECONMY, 

IOUT, 

ICOL 



17 

LOGICAL ECONMY 





18 

DOUBLE PRECISION 

RWT , TOTWT 

.WEIGHT 



19 

DATA INTER/3/, INPUT5/5/ , 

LIST/6/ 



20 

DATA ( FMTTRI I I ) , I 

=1,A)/6H(5H RO, 6HW I 

5,2, 6HX , I 8G 1 , 

6H5.6I) / 

21 


COMMON/MAX/M AX PLT 22 
MAXPLT SHOULD BE THE NUMBER OF SINGLE LENGTH WORDS IN COMMON/BIG/ 23 
BEGINNING AT THE FIRST LOCATION OF SUMXY 2A 
DATA MAXPLT/10700/ 25 
END 26 


BFTC NEWRAP 


THIS IS NEWRAP, MAIN PROGRAM FOR REGRESSION ANALYSIS PROVIDING 
INTERNAL EVALUATION OF RESULTS. 

COMMON / FRMTS/ FMT( 13), FMTTRI ( 14) 

COMMON /BIG/BI 60,9) , SUMXY ( 60, 9) ,SUMXX< 1830) , SUMXX I ( 1 830 ) 

DOUBLE PRECISION B, SUMXY .SUMXX, SUMXX I 

COMMON /MED/BOI 9) .SUMXI69 ), SUMX2169I , SUMY2I9) , ZEANI69) , 

X CONI 99) .ERRMSI9) , I DENT 1 13) , I DOUT ( 60 ) , NCON I 200) , NTERMI 69) , 

X NTRAN ( 100 ) , NXCOD (100), POOL ED I 9 ) , R EP VAR I 9 ) , RE SMS t 9 ) , X I 99 ) 
DOUBLE PRECISION BO.SUMX.SU MX2 , SUMY2 , Z E AN 
COMMON/SMALL/ BYPASS, BZERO 

X IFTT , IFWT, 

X ISTRAT, JCOL » 

X NERROR, NODEP, 

X NOVAR, NP D EG, 

X P, PREDCT , 

X STORY I , STORYC, 

X ERRFXD, ECONMY , IOUT, ICOL 
LOGICAL ECONMY 

DOUBLE PRECISION RWT , TOT WT , W E I GHT 


LOGICAL 
XIFSSR, 
XSTORYC . 
LOGICAL 
LOGICAL 


BYPASS, 

IFTT, 

STORYX, 

XSAVE 

PNCH 


BZERO, 

IFWT, 


DIMENSION XCHKI60) 


I-************************************* ******************** ************* 


EQUIVALENCE ( NAR AY, SUMMX I ) , ( S , BO ) , l S SQ , REP VAR ) 
DIMENSION NARAYI 1830 ) , S ( 9 ) , S SQ I 9 ) 


ZERO OUT ALL DATA ARRAYS EACH NEW DATA SET 
100 DO 101 J=1 ,4740 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 


1, DELETE, FIRST, IFCHI, 

IFSSR, 

13 

INPUT, 

INPUT5, 

INTER, 

1A 

KONNO, 

LENGTH, 

LIST, 

15 

NOOB, 


NOTERM, 

16 

NRES, 

NTRANS, 

NWHERE, 

17 

REPS, 

RWT, 


18 

STORYX, 

TOTWT, 

WEIGHT, 

19 




20 




21 

GHT 



22 

DELETE 

, IFCHI, 


23 

REPS, 

PREDCT 

* 

2 A 

, FIRST 

.ERRFXD 


25 


26 

27 

28 

29 

30 

31 

32 

33 
3 A 

35 

36 
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101 B ( J 1 1 ) =0 . 0D0 37 

DO 102 J=l,225 38 

102 80( J ) = 0. ODD 39 

C 40 

C***************************** ****************************************** 41 

C READ IDENTIFICATION CARD AND OPTIONS CARD 42 

C 43 

READ(INPUT5,110) I » I DENT 44 

WRITE! LIST, 111) IDENT 45 

FIRST=.TRUE. 46 

ERRFXD=. FALSE. 47 

113 I F { I > 120,120,115 48 

115 R£A0(INPUT5,300) FMT 49 

WRITE! LIST, 301) FMT 50 

1=1-1 51 

GO TO 113 52 

120 READ(I NPUT5, 1282 ) NOV AR, NODEP, NOT ERM, NOOB, NTKEEP 53 

WRITE! LIST, 1283) NOV AR, NODEP .NOTERM, NOOB 54 

IF!NT<EEP.NE.O) WR I TE { L I ST , 1 307 ) NTKEEP 55 

RE AD (I NPUT5.117) BZERO , IFTT , IFWT, IFCHI, STORYX, IFSSR, 56 

X ECONM Y , ISTRAT, PNCH 57 

WRI TE! LIST, 118) BZERO, IFTT, IFWT, IFCHI, STORYX, IFSSR, ECONMY , 58 

XISTRAT, PNCH 59 

C 60 

C*********************************************************************** 61 

C THESE ARE INITIALIZATIONS MADE BEFORE EACH SET OF DATA 62 

C ICOL DETERMINES THE NUMBER OF VARIABLES READ PER OBSERVATION 63 

C JCOL IS THE NUMBER OF TERMS IN THE TOTAL REGRESSION EQUATION 64 

C LENGTH IS THE NUMBER OF STORES NEEDED FOR THE MATRICES 65 

LENGTH= NOTE RM* I NOTERM+1 )/2 66 

ICOL=NOVAR + NODEP 67 

JCOL = NOTERM +NODEP 68 

NWHERE= NOTERM 69 

REWIND INTER 70 

DO 143 J= 1 , 60 71 

I DOUTI J ) = J 72 

140 NTE RM! J ) = J 73 

DO 145 J=l,100 74 

NXCOD! J ) = J 75 

NTRAN! J ) =0 76 

145 NCON ( 2 *J ) = J 77 

C 78 

C*********** ********** ******** ****************************************** 79 

IF 1 BZERO) WR I TE ! L 1ST , 190 ) 80 

IF!. NOT. BZERO) WR I TE ! L I ST, 170 ) 81 

C**** ******************************************* ************************ 82 

READ! I NPUT5.282) NTR ANS, KONNO 83 

IF(NTRANS.EQ.O) GO TO 255 84 

220 READ IINPUT5,230)(NTERM!K),K=1,JC0L) 85 

WRITE! LIST, 235) 1 NTERMIK ) , K= 1 , JCOL ) 86 

READ ( INPUT5,230)(NXCOD(I),NTRAN!I),NCON(2*I-1),NCON!2*I) , 1=1, NTR 87 

AANS ) 88 

WRITE! LIST, 240) ! NXCOD! I ), NTRAN! I ) , NCON ! 2* 1-1 ) ,NC0N(2*I ) ,1=1, 89 

X NTR ANS) 90 

IF ( KONNO) 255,255,250 91 

250 READ < INPUT5, 260) ICON! I ), 1 = 1, KONNO) 92 

WRITE! LIST, 262) I CON ( I ) , 1= 1 , KONNO ) 93 

C**** ************ ********* ••••••*•*••*•••**••••*••••••**••••••••*• »»»»•• 94 

C 95 

255 READ!INPUT5,257) DELETE, P 96 

I F { DELETE ) I FTT=. TRUE . 97 

C 98 
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noon 


* * * * 


IF THERE ARE REPLICATED POINTS READ IN THE NUMBER OF POINTS AND 
THE NUMBER OF REPLICATIONS. SINGLE DATA POINTS ARE DATA POINTS 
REPLICATED ONCE. OBSERVED DATA MUST BE ARRANGED IN THE ORDER 
IMPLIED HERE. 

READ! I NPUT5, 257 ) REPS 
XSAVE=. FALSE. 

265 IF{. NOT. REPS) GO TO 290 

RE AD ( I NPUT5.282 ) I RE P , ( N AR AY (I) , I = 1, I REP ) 

WRITE! LIST, 284) IREP 

WRI TE( LIST, 233) (NARAY! I ) ,1 = 1, IREP) 

NPDC G= 0 

I RE P = 1 

IC = NAR AY! 1 ) 

XSAVE= • TRUE. 

DO 315 1=1* NODEP 
POOLED! I )=0.0 
S t I ) =3. 0 

315 SSO(I)=0.0 


C 


C READ VARIABLE FORMAT FOR DATA 

290 READ! I NPUT5, 110) INPUT, FMT 
WRITE(LIST,lli) FMT 
310 TOT WT = 0. ODO 
WEIGHT = 1.000 
WR I TE t L I ST, 301 ) IDENT 


C 


C READ IN INPUT VARIABLES 

DO 493 J=1,nOOB 
330 IF! .NOT. IFwT 1 GO TO 350 
340 READ (INPUT, FMT) ( X ( I ) , I =1 , I COL ) 

GO TO 360 

350 READ ( INPUT, FMT) ( X ( I ) ,1=1, ICOL) , WEIGHT 
360 CONTINUE 

I F ( ECONMY ) WRITEILIST ,381) J , ( X ( I ) , I = 1 , I COL ) 
381 FORMAT ( 1H I 4 , 9G 1 4. 6/ ( 5X , 9G 14 . 6 ) ) 
IF(NTRANS.EQ.O) GO TO 450 
IFIECONMY) GO TO 390 
WRITE! LI ST, 3 70) WEIGH T,J 
WRITE (LIST, 380) (X(I ), 1=1, ICOL) 

390 CALL TRANS 
420 DO ^33 K= 1 , JCOL 

I=NTERM(K) 

X ( K ) = CON! I ) 

430 CONTINUE 
450 CONTINJE 

IF ( ECDNMY ) GO TO 4609 
WRITE! LIST, 460) J 

461 WRITE (LIST, 380) (X! I ), 1 = 1, JCOL) 

4609 CONTINJE 

IF(IFOHI) WRITE! INTER) ( X ( I ) , I = 1 , 69 ) , WE I GH T 
IF! .NOT.XSAVE) GO TO 4611 
DO 4610 K = l, NOTERM 

4610 XCHK { < ) =X ( K ) 

XSAVE=. FALSE. 

4611 CONTINUE 


C 


i » * * » *■ * i 


l»*»**#******#**»»*»***»»**s*** 


99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 
1 54 

155 

156 

157 

158 
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C COMPJTE THE ERROR VARIANCE FROM REPLICATED DATA 159 

IF(. NOT. REPS) GO TO 480 160 

IGOTO =1 161 

IF(NARAY( IREP1.GT.1) IGOTO=2 162 

IF(J.SE.IC) WRITE ( 6, 462 ) IREP 163 

DO 475 1 = 1 » NODEP 164 

IF(I-l) 4629,4629,464 165 

4629 DO 463 K=l, NOTERM 166 

IF(X(O.NE.XCHK(K) ) GO TO 2001 167 

463 CONTINUE 168 

464 CONTIVUE 169 

KBAR=VOTERM+I 170 

S(I)=S( D+XIKBAR) 171 

SSQ( I) =SSQ( I )+X(KBAR )**2 172 

IF(J-IC) 475,465,465 173 

465 GO TO (468,466) , IGOTO 174 

466 ZEAN(l)=S(I)/FLOAT(NARAY(IREP) ) 175 

SSQ(I) = SSQ(I) - ZE AN ( 1 ) *S ( I ) 176 

POOLED (I)=POOLED(I)+SSQ( I) 177 

WRITE( LIST, 467) I , SSQ ( I ) ,S( I ) ,ZEAN( 1 ) 178 

468 IF t I. LT. NODEP) GO TO 469 ■ 179 

NPDEG=NPDEG+NARAY( IREP ) -1 180 

I REP=I REP+1 181 

IC = IC + NARAY ( IREP ) 182 

WRITE(LIST,4671) 183 

469 S( I ) =3.0 184 

SSQ ( I ) =0, 0 185 

XSAVE= • TRUE, 186 

475 CONTIVUE 187 

C 188 

C*********************************************************************** 189 

C CALCULATE SUMS, SUMS OF SQUARES AND SUMS OF CROSS PRODUCTS. 190 

480 CALL SUMUP 191 

490 CONTIVUE 192 

C 490 CONTINUE IS THE END OF THE LOOP FOR READING DATA CARDS 193 

IF(. NOT. REPS) GO TO 496 194 

DO 493 1=1, NODEP 195 

REPVAR(I)=POOLEO( I ) / FLOAT ( NP DEG ) 196 

493 CONTIVUE 197 

496 CONTIVJE 198 

C 199 

C*********************************************************************** 200 

C ALL DATA HAS BEEN READ IN AND THE XTRANSPOSEX AND XTRANSPOSEY 201 

C MATRIX HAVE BEEN CALCULATED. 202 

C NOW WRITE THE MATRICES 203 

CALL MFIX 204 

REWIND INTER 205 

GO TO 640 206 

C 207 

C*************** ************************************************* ******* 208 

C THIS CODING DELETES THE DATA FROM THE SUMXX MATRIX 209 

C CORRESPONDING TO THE INDEPENDENT TERM DELETED 210 

6500 CONTIVUE 211 

I R= I OJT-1 212 

IC= NDTERM - I OUT 213 

IF ( IC.EQ. 0) 50 TO 6700 214 

I NOCH= I OUT * I R/2 215 

INEW = INOCH 216 

IOLD a INEW + I OUT 217 

IRC=0 218 

IBC=0 219 

ITC=0 220 
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DO 6630 I=I0L0, LENGTH 221 

INEW = INEW+1 222 

IOLD=I OLD + 1 223 

IF! ITC.GT.O) GO TO 6540 224 

I RC= 1 3 C + 1 225 

IF! IRS.GT.IR) GO TO 6530 226 

SUMXX! INEW)=SUMXX! IOLD) 227 

GO TO 6600 228 

6530 I BC= I B C + 1 229 

ITC = IBC 230 

IOLD = IOLD+1 231 

I RC= 3 232 

6540 ITC = ITC -1 233 

SUMXX! INEW) = SUMXX( IOLD) 234 

6600 CONTI 'IDE 235 

6700 LENGTH = LENGTH-NOTERM 236 

N0TERM= NOTERM -1 237 

JC0L= NOTERM+NODEP 238 

239 

240 

INVERT THE SUMXX MATRIX AND COMPUTE REGRESSION COEFS 241 

AND S J MS OF SQUARES DUE TO REGRESSION IN THE MATRIX INVERSION 242 

ROUTINE 243 

640 CONTINUE 244 

CALL MATINV 245 

FIRST*. FALSE. 246 

247 

************************************************************************ 248 

WRITE(XTX) INVERSE. THIS MATRIX TIMES ERROR MEAN SQUARE (ERRMS) 249 

IS THE VARIANCE-COVARIANCE MATRIX OF REGRESSION COEFFICIENTS. 250 

IF(FCONMY) GO TO 970 251 

WRITE! LIST, 700) 252 

CALL TRIANGU, SUMXXl, NOTERM, 8, FMTTRI, 2) 253 

254 

255 

IF A VARIABLE HAS BEEN DELETED ADJUST COUNTERS AND RECOMPUTE THE 256 

REGRESSION. IF NO VARIABLE HAS BEEN DELETED CONTROL WILL PASS 257 

FROM THE TTEST ROUTINE TO THE CHI-SQUARE OPTION. 258 

970 CONTINUE 259 

I F ( . NOT . I FTT ) GO TO 1020 260 

980 WRITE (LIST, 301) IDENT 261 

CALL TTEST($1020,NTKEEP) 262 

I F ( NOOEP— 1 ) 985,990,985 263 

985 WRITE! LIST, 986) NODEP 264 

NODEP= 1 265 

990 J= JCOL -1 266 

DO 995 K= I OUT , J 267 

NTERM! K)=NTERM(K+1) 268 

ZEANIO = ZE AN ( K+l ) 269 

SUMX(K) = SUMX ( K+l ) 270 

SUMX2 ( K ) = SUMX2 ( K+1-) 271 

IDOUT(K) = I DOUT ( K+l ) 272 

SUMXY! K , 1 ) = SUMX Y ( X+ 1 , 1 ) 273 

995 CONTINUE 274 

IF(NOTERM.EQ.l) GO TO 1000 275 

GO TO 6500 276 

1000 WRITE! LIST, 1005) 277 

NOTERM=0 278 

GO TO 1035 279 

280 

*********************************************************************** 281 

282 
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on on 


1020 I F ( .NOT. IFCHI ) GO TO 1035 283 

1030 WRITE( LIST, 301) I DENT 284 

CALL 3UTPLT { PNCH ) 285 

286 

287 

1035 READ( INPUTS, 117) PREDCT 288 

I F ( . NOT. PREDCT ) GO TO 100 289 

CALL PREDIC 290 

1040 GO TO 100 291 

292 

*********** 293 

2001 WRITE( LIST, 1306) 294 

STOP 295 

C *************************** ******************************************** 296 

8001 FORMAT ( 1H1 ) 297 

8002 FORMAT ( 1H2 ) 298 

110 FORMAT (12, 13A6 ) 299 

111 FORMAT ( 1 Hi , 13A6 , A2 ) 300 

117 FORMAT (7L1, II, LI ) 301 

118 FORMAT ( 1H 7L1,I1,L1) 302 

170 FORMAT ( 33H THERE IS NO BO TERM IN THE MODEL) 303 

190 FORMAT ( 26H THERE IS A BO TO ESTIMATE ) 304 

230 FORMAT (4012) 305 

235 FORMAT ( 11H NTERM(K)= /(1H 3014)) 306 

240 FORMAT (25H THE TRANSFORMATIONS ARE /(1H 5(414, 5X) ) ) 307 

257 FORMATt 1L1, F3.3) 308 

260 FORMAT (5E15.7) 309 

262 FORMAT ( 19H TH& CONSTANTS ARE /((1H 8G15.7))) 310 

282 FORMAT ( 2014 ) 311 

283 FORMAT ( 1H 2014) 312 

284 FORMAT ( 11H THERE ARE I5.16H REPLICATE SETS ) 313 

300 FORMAT ( 13A6, 1A2 ) 314 

301 FORMAT ( 1H 13A6.A2) 315 

370 FORMAT (1H0»29H0BSERVED VARIABLES, WEIGHT = G14.6, 6X, 15H0BSERVATI ON 316 

1 = ,15 ) 317 

380 FORMAT ( 1H 9G14.6) 318 

460 FORMAT ( 1H , 37HTERMS OF THE EQUATION, OBSERVATION = ,15) 319 

462 FORMAT ( 18HK** REPLICATE SET 1 5, 3X» 100( 1H» ) ) 320 

4671 FORMAT ( 1H 125(1H*)) 321 

467 FORMAT ( 14H DEP. VAR. 16, 8H SSQ=G14.7,8H SUM=G14.7,8H M 322 

X EAN= 314.7) 323 

540 FORMAT ( 1H 8G14.7) 324 

560 FORMAT ( 21H2X TRANSPOSE X MATRIX ) 325 

670 FORMAT ( 25 H2 COR RELATION COEFFICIENTS ) 326 

700 FORMAT ( 32H2 ( X TRANSPOSE X) INVERSE MATRIX ) 327 

986 FORMAT ( 39H THE NUMBER OF DEPENDENT VARIABLES WAS I3,83H IT IS BE 328 

XING SET TO ONE AND THE REJECTION OPTION EXERCISED ON DEPENDENT VAR 329 

XIABLE 1 ) 330 

1005 FORMAT ( 39H THERE IS NO EVIDENCE OF A REGRESSION. / 331 

X 74H USE THE MEAN RESPONSE FOR THE BEST ESTIMATE OF THE DEPEND 332 

XENT VARIABLE (S) • ) 333 

1282 FORMAT (314,15,14) 334 

1283 FORMAT ( 1H 314,15) 335 

1306 FORMAT ( 40H REPLICATE SETS ARE NOT GROUPED PROPERLY ) 336 

1307 FORMAT ( 11H THE FIRST I2.64H TERMS OF THE MODEL WILL BE RETAINED R 337 

XEGARDL ESS OF SIGNIFICANCE ) 338 

END 339 
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on on 


SIBFTC MATINX 


C 
C 
C 

C*** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


SUBROJTI NE MATINV 


•if****************************************************************** 


(X-TRANSPOSE X) 


PURPOSE 

1) COMPUTE EIGENVALUES AND EIGENVECTORS OF 

MATRIX IF REQUESTED. ( STORYX= . TRUE . ) 

2) COMPUTE (X TRANSPOSE X) INVERSE 

3) COMPUTE REGRESSION COEFFICIENTS 

A) COMPUTE OTHER REGRESSION STATISTICS 

SUBROUTINES NEEDED 
BORD 
LOC 
EIGEN 
INVXTX 
RECT 
RSTATS 
TRI ANG 


REMAR<S 

THE EIGENVALUES ARE COMPUTED AS AN AID IN DETERMINING THE 
CONDITION OF THE SYSTEM OF EQUATIONS FOR THE REGRESSION 
COEFFICIENTS. EXAMINATION OF THEM AND THEIR ASSOCIATED 
EIGENVECTORS MAY SHOW THAT CERTAIN SETS OF INDEPENDENT 
VARIABLES ARE HIGHLY CORRELATED AND NOT EASILY LIABLE TO 
INDEPENDENT STUDY. 

SUBROUTINE MATINV 


in,**************************************************************#****** 

C0MM0N/BIG/B(60,9) ,SUMXY (60,9) , SUMXXI 1830) .SUMXXI ( 1830) 

DOUBLE PRECISION B, SUMXY , SUMXX, SUMXX I 

COMMON /MED/ BO( 9) ,SUMX< 69 ),SUMX2<69) ,SUMY2( 9) , ZEAN(69> , 

X CONI 99) ,ERRMS(9) ,1 DENT (13) , I DOUT { 60 1 , NCON ( 200 ) ,N TERM (69) , 

X N IRANI 100) , NXCOD (100), POOL ED ( 9 ) , REP VAR ( 9 ) , RE SMS ( 9 ( , X ( 99 ) 

DOUBLE PRECISION BO,SUMX,SU MX2 , SUM Y2 , ZE AN 


XSTORYC. STORYX, STORYI, 

DOUBLE PRECISION RWT . TOT WT .WEIGHT 
DIMENSION AI1),C(1),XTX(3) 

EQUIVALENCE (A, SUMXX I ) , ( C, SUMXX I C 915 ) ) 
DOUBLE PRECISION SUM 

DATA (XTX( I ) , 1 = 1,3) / 6HX TRAN, 6HSP0SE 


FIRST , ERRFXD 


6HX 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 
11 
12 
13 
1 A 

15 

16 

17 

18 

19 

20 
21 
22 
23 

2 A 

25 

26 

27 

28 

29 

30 

31 

32 

33 

3 A 

35 

36 

37 

38 


COMMON /FRMTS/ 

FMT ( 13) , FMTTRI ( 1A) 
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COMMON/SMALL/ 

BYPASS, BZERO 

.DELETE, FIRST, IFCHI, 

IFSSR, 

AO 

X I FTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

A1 

X ISTRAT, 

JCOL , 

KONNO, 

LENGTH, 

LIST, 

A2 

X NERROR, 

NOOEP, 

NOOB, 


NOTERM, 

A3 

X NOVAR, 

NPDEG, 

NRES, 

NTRANS , 

NWHERE , 

AA 

X P, 

PREDCT , 

REPS, 

RWT, 


A5 

X STORYI, 

STORYC, 

STORYX, 

TOTWT, 

WEIGHT, 

A6 

X ERRFXD, ECONMY 

, I OUT , ICOL 




A7 

LOGICAL ECONMY 





A8 

LOGICAL BYPASS, BZERO, 

DELETE 

, IFCHI, 


A9 

XIFSSR, I FTT 

, IFWT, 

REPS, 

PREDCT 

, 

50 


*********************************************************************** 

IORDER = NOTERM 

IF ( NOTERM-1 ) 10,10,12 


51 

52 

53 
5 A 

55 

56 

57 

58 

59 
A n 


53 


onno non non oo 


10 SUMXX I ( 1 ) = 1 .U/SUMX2 ( 1 ) 61 

GO TO 350 62 

63 

64 

12 I F ( .NOT • STORYX) GO TO 30 65 

DO 14 1=1, LENGTH 66 

A ( I ) =S UMXX ( I ) 67 

14 CONTINUE 68 

16 CALL E IGEN! A,C ,I0RDER,0) 69 

WRITE! LIST, 17) (XTX(I) , 1=1,3) 70 

J=0 71 

DO 18 1=1 , IORDER 72 

J=J+I 73 

18 A! I ) =A ( J ) 74 

WRI TE( L 1ST , 19 ) ( A ( I ) , I = 1 , IORDER ) 75 

WRITE! LIST, 20) 76 

CALL RECTI IORDER, IORDER, IORDER, IORDER, C,X »FMTTR 1,1) 77 

30 DO 35 1=1, LENGTH 78 

35 SUMXXI ( I )=SUMXX! I ) 79 

80 

*********************************************************************** 81 

NO S JBMODELS TO ANALYZE SO INVERT A DIRECTLY BY GAUSS 82 

49 IF! IFSSR) GO TO 50 83 

CALL I NVXTX ( SUMXX I , NOTERM, D, 1.0 ) 84 

GO TO 60 85 

86 

*********************************************************************** 87 

C SUBMODELS HAVE BEEN REQUESTED SO WE USE BORDERING 88 

50 IOROER=0 89 

55 IOROER=IORDER +1 90 

CALL BORD! IORDER, SUMXXI ) 91 

60 CONTINUE 92 

93 

*********************************************************************** 94 

COMPUTE COEFFICIENTS AND PRINT THEM 95 

96 

350 DO 373 J=1,N00EP 97 

DO 373 K=1 , IORDER 98 

B!K,J)=O.ODO 99 

DO 373 L=l, IORDER 100 

CALL LOC(L,K»IR) 101 

B(K,J> = B(K,J) * SUMXXI ( I R ) *SUMX Y { L , J ) 102 

370 CONTINUE 103 

C 104 

WRITE! LIST, 380) IDENT 105 

WRI TE(LIST,382) 106 

IF! .N3T.BZER0) GO TO 400 107 

DO 393 J= 1 , NODEP 108 

SUM=0. ODO 109 

KBAR = NOTERM + J 110 

DO 385 K=l, IORDER 111 

SUM = SUM + B { K, J ) *Z E AN! K ) 112 

385 CONTINUE 113 

BO ! J ) = ZEAN(KBAR) -SUM 114 

390 CONTINUE 115 

WRITE(LIST,395) 116 

WRITE! LIST, 397) { BO ! K ) , K=1 , NODE P ) 117 

400 WRITE! LIST, 410 ) 118 

DO 433 J=l, IORDER 119 

WRITE! LIST, 432) I DOUT ! J ) , l B ( J , K ) , K= 1 , NODEP ) 120 

430 CONTINUE 121 

C 122 
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123 

COMPUTE REGRESSION STATISTICS IN R STATS 124 

125 

CALL RSTATSI IORDER) 126 

127 

*********************************************************************** 128 

IF IORDER IS LESS THAN NOTERM WE HAVE USED THE BORDERING OPTION 129 

AND MUST GO BACK TO FINISH. 130 

131 

I F C IORDER-NGTERM) 55,500,500 132 

500 STORYK=. FALSE. 133 

IFSSR=. FALSE. 134- 

RETURN 135 

17 FORMAT (34H2THE FOLLOWING ARE EIGENVALUES OF 2A6.A1, 7H MATRIX) 136 

19 FORMAT ( 1H 8G16i7) 137 

20 FORMAT ( 132H2THIS IS THE MODAL MATRIX OR MATRIX OF EIGENVECTORS. El 138 

1GENVEC TORS ARE WRITTEN IN COLUMNS LEFT TO RIGHT IN SAME ORDER AS 139 

2EIGENVALUES ) 140 

380 FORMAT (1H1,13A6,1A2) 141 

382 FORMAT! 61H EACH COLUMN CONTAINS THE COEFFICIENTS FOR ONE DEPENDEN 142 

XT T C RM ) 143 

395 FORMAT ( 20H CONSTANT TERM (BO) ) 144 

397 FORMAT (4X.9G14. 6) 145 

410 FORMAT ( 36H REGRESSION COEFFICIENTS (B1,...,BK) ) 146 

432 FORMAT ( 1H I3.9G14.6) 147 

END 148 


IBFTC TTt-STX 


*********************************************************************** 

SUBROUTINE TTEST 
PURPOSE 

COMPUTE THE T-STATISTICS FOR EACH REGRESSION TERM AND 
ITS TWO TAILED SIGNIFICANCE LEVEL. THEN DETERMINE THE 
TERM WITH THE LEAST SIGNIFICANCE AND RETURN THIS 
INFORMATION TO NEWRAP 

*••••*••••»•••*»**•»••••••••*•••*•••»•••••»•••*•••••*•••*•••••••••••••• 

SUBROUTINE TTEST (*, NTKEEP ) 

***************************** ******** ********************************** 

COMMON / FRMTS/ FMT I 1 3 ) , FMTTR I I 14 ) 

COMMON ZB1 G/B (60,9 ) , SUMXY (60,9) ,SUMXX( 1830) , SUMXXI ( 1830) 

X , DUMMY ( 1 ) 

DOUBLE PRECISION B, SUMXY , 5UMXX , SUMXX I 

CQMM0N/MED/BQ(9),SUMX(69),SUMX2(69)» SUMY 2(9) » Z t AN ( 69 ) , 

X CON( 99) , ERRMS ( 9 ) , I DENT ( 13) , I DOUT ( 60 ) , NCON ( 200 ) , NTERM ( 69 ) , 

X NTRAN( 100 ) ,NXC0D(100) , POOLED! 9 1 , REP VAR l 9 ) , Rt SMS I 9 ) ,X(99) 

DOUBLE PRECISION BO,SUMX,SU MX2 , SUM Y 2 , Z E AN 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 


COMMON/SMALL/ 

BYPASS 

,BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

23 

X 

I FTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

24 

X 

ISTRAT, 

JCOL, 

KONNO, 

LENGTH, 

LIST, 

25 

X 

NERROR, 

NOUEP, 

NOOB , 


NOTERM, 

26 

X 

NOVAR, 

NP DEG, 

NRES , 

NTRANS , 

NWHERE , 

27 

X 

P, 

PREUCT 

, REPS, 

RWT , 


28 
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STORYC, STORYX, 

I OUT , I COL 


X SIORYI, 

X ERRFXD, ErCONMY, 

LOGICAL ECONMY 
LOGICAL SATRTD 

LOGICAL BYPASS. BZERO, 

XIFSSR, IFTT, IF WT . 

XSTORYC, STORYX, STORYI, 

DOUBLE PRECISION RWT , TOT WT , W E I GHT 


TOTWT, 


WEIGHT, 


DELETE, IFCHI, 

REPS, PREDCT , 

FIRST , ERRFXD 


LOGICAL MAKENU, NOZERO 
DIMENSION T(35,13) .PLEVELI 13) 

DIMENSION DEVB ( 60,9) , PR0B(60,9) ,TT(60,9) 

EQUIVALENCE ( DUMM Y ( 9 1 ) , D EVB , TT ) , I DUMMY ( 650 ) , PROB ) 
EQUIVALENCE (P.PWANT) 


DATA (PLEVELI JJ ) ,JJ=1, 13) / 0. 10, 0. 20, 0 . 30 , 0. 40, 0. 50 , 0, 60, 0. 70 , 
10.80,3.90,0.95,0.98, 0.99,0.999 / 

DATA (T (1,JJ) , JJ= 1 » 13) /O. 158, 0.325, 0.510, 0.727, 1.000, 1.376, 

1 1.963,3.078,6. 314, 12.706,31.821,63.657,636.619 /, 

2 ( T ( 2 , J J ) »JJ=1, 13 ) /O. 142, 0.289, 0.445, 0.617, 0.816, 1.061, 

3 1.386,1.886,2.920,4.3027,6.965,9.925,31.598 /, 

4 (T(3»JJ),JJ=1,13) /O. 137, 0.277, 0.424, 0.584, 0.765, 0.978, 

5 1.250,1.638,2.353,3.1825,4.541,5.841,12.924 /, 

6 IT (4,JJ),JJ=1,13) /O. 134,0. 271, 0.414, 0.569, 0.741, 0.941, 

7 1.190,1.533,2.132,2.7764,3.747,4.604,8.610 /, 

8 (T(5,JJ),JJ=1,13) /O. 132, 0.267, 0.408, 0.559, 0.727, 0.920, 

9 1.156,1.476,2.015,2.5706,3.365,4.032,6.869 /, 

A (T (6,JJ),JJ=1,13) /O. 131, 0.265, 0.404, 0.553, 0.718, 0.906, 

B 1.134,1.440,1.943,2.4469,3.143,3.707,5.959 /, 

C (T(7, JJ) , JJ=1 , 13 ) /O. 130, 0.263, 0.402, 0.549, 0.711, 0.896, 

D 1.119,1.415,1.895,2.3646,2.998,3.499,5.408 /, 

E ( T(8, JJ ) , JJ=1 ,13) /O. 130, 0.262, 0.399, 0.546, 0.706, 0.889, 

F 1.108,1.397,1.860,2.3060,2.896,3.355,5.041 /, 

G (TI9.JJ) ,JJ=1»13) /O. 129, 0.261, 0.398, 0.543, 0.703, 0.883, 

H 1.100,1.383,1.833,2.2622,2.821,3.250,4.781 /, 

I IT (10, JJ) ,JJ=1,13) /O. 129,0. 260, G. 397, 0.542, 0.700, 0.879, 

J 1.093,1.372,1.812,2.2281,2.764,3.169,4.587 / 

DATA (T(11»JJ),JJ=1,13) /O. 129, 0.260, 0.396, 0.540, 0.697, 0.876, 

1 1.088,1.363,1.796,2.2010,2.718,3.106,4.437 /, 

2 (T ( 12, JJ) »JJ = 1»13) /0. 128, 0.259, 0.395, 0.539, 0.695, 0.873, 

3 1.083,1.356,1.782,2.1788,2.681,3.055,4.318 /, 

4 (T( 13, JJ) , JJ= 1 , 13 ) /O. 128, 0.259, 0.394, 0.538, 0.694, 0.870, 

5 1.079,1.350,1.771,2.1604,2.650,3.012,4.221 /, 

6 (T(14,JJ),JJ=1,13) /O. 128, 0.258, 0.393, 0.537, 0.692, 0.868, 

7 1.076, 1.345,1.761, 2.1448,2.624,2.977,4.140 /, 

8 (T(15,JJ)»JJ=1,13) /O. 128, 0.258, 0.393, 0.536, 0.691, 0.866, 

9 1.074,1.341,1.753,2.1315,2.602,2.947,4.073 /, 

A (T(16,JJ),JJ=1,13) /O. 128, 0.258, 0.392, 0.535, 0.690, 0.865, 

B 1.071,1.377,1.746,2.1199,2.583,2.921,4.015 /, 

C (T(17»JJ) ,JJ=1,13» /O. 128,0.257,0.392,0.534,0.689,0.863, 

D 1.069,1.333,1.740,2.1098,2.567,2.898,3.965 /, 

E (T(18,JJ),JJ=1,13) /O. 127, 0.257, 0.392, 0.534, 0.688, 0.862, 

F 1.067,1.330,1.734,2.1009,2.552,2.878,3.922 /, 

G (T ( 19, JJ) »JJ=1,13 ) /O. 127, 0.257, 0.391, 0.533, 0.688, 0.861, 

H 1.066,1.328,1.729,2.0930,2.539,2.861,3.883 /, 

I (T (20, JJ ) ,JJ=1,13) /O. 127, 0.257, 0.391, 0.533, 0.687, 0.860, 

J 1.064,1.325,1.725,2.0860,2.528,2.845,3.850 / 

DATA (T(21,JJ), JJ=1, 13) / 0 . 12 7, 0. 25 7, 0 . 391 , 0 . 532 , 0 . 686 , 0. 8 5 9 , 

1 1.063,1.323,1.721,2.0796,2.518,2.831,3.819 /, 

2 (T(22,JJ),JJ=1,13) /O. 127,0. 256, C.J90, 0.532, 0.686,0. 858, 


29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 
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3 

4 

5 

6 

7 

8 
9 
A 
B 
C 
D 
E 
F 
G 
H 

1 
J 

DATA 

1 

2 

3 

4 

5 

6 
7 


1.061, 1.321,1.717, 2.0739,2.508,2.819,3.792 /, 

(T (23, JJ) ,JJ=1, 13) /0. 12 7, 0.256, 0.390, 0.532, 0.685,0. 858, 
1.06G, 1.319,1. 714,2.0687,2.500,2.807,3.767 /, 

(T(24,JJ),JJ=1,13) /0. 127, 0.2 56,0.390, 0.53 1,0. 685,0. 857 , 
1.059,1.318,1.711,2.0639,2.492,2.797,3. 745 /, 

{T< 25, JJ) , JJ=1, 13 ) / 0.1 2 7, 0.256, 0.390, 0.53 1,0.684,0. 856, 
1.058,1.316,1.708,2.0595,2.485,2.787,3.725 /, 

(T(26,JJ),JJ=1,13> /0. 127, 0.2 56, 0.3 90, 0.531, 0.684, 0.856, 
1.058, 1.315,1. 706, 2. 05 55 , 2. 479, 2 . 779 , 3. 707 /, 

(T (27»JJ),JJ=1»13) /0. 12 7, 0.2 56, 0.389, 0.531, 0.684, 0.85 5, 


1.057, 
( T ( 28 
1.056, 
( T ( 29 
1.055, 
(TOO 
1.055, 
(T C 31 
1.050, 
( T ( 32 
1.046, 
( T C 33 
1.041, 
( T ( 34 
1.036, 


1.314, 1 
, J J ) , J J 

1.313.1 
, JJ) , JJ 

1.311.1 
, JJ ) , JJ 

1.310.1 
,JJ) , JJ 

1.303.1 
, JJ I , JJ 

1.296.1 
,JJ) , JJ 

1.289.1 
, J J ) , J J 

1.282.1 


. 703,2. 
= 1,13 ) 

. 701,2. 
=1,13) 

, 699,2. 
=1,13) 
,697, 2. 
=1,13) 
,684, 2. 
=1,13) 
.671, 2. 
=1,13) 
.658, 1. 
=1,13) 

. 645. 1 . 


0518. 2. 
70.127, 

0484.2. 
/ 0 . 1 27, 

0452. 2. 
/ 0 . 1 27 , 

0423.2. 
/ 0 • 1 26, 

0211 . 2 . 

/0.126, 

0003.2. 
/O. 126, 

9799. 2. 
/0.126, 

9600.2. 


tlJli. I I 1 » .3 • ( 


0.256, 

467.2. 
0.256, 

462.2. 
0.256, 

457.2. 
0.255, 

423.2. 
0.254, 

390.2. 
0.254, 

358.2. 
0.253, 

326.2. 


0.389, 

763.3. 
0.389, 

756.3. 
0.389, 

750.3. 
0.388, 

704.3. 
0.387, 

660. 3. 
0.386, 

617.3. 
0.385, 

576. 3. 


674 

0.530,0. 

659 

0. 530 , 0. 
646 

0.529,0. 

551 

0.527,0. 

460 

0.526,0. 

373 

0.524,0. 

291 


7 , 

683.0. 855, 
7, 

683.0. 854, 
7, 

683.0. 854, 
7 

681.0. 851, 
/ * 

679.0. 848, 
7, 

677.0. 845, 

/ , 

674.0. 842, 
7 


T(II,JJ) IS THE T-STATISTIC AT THE TABULATED DEGREES OF FREEDOM 
(II) AND AT THE TABULATED PROBABILITY LEVELS (JJ). 

It=3EGREFS OF FREEDOM, EXCEPT FOR 

11=31 IS FOR 40 DEGREES 

11=32 IS FOR 60 

11=33 IS FOR 120 

11=34 IS FOR INFINITY 


JJ 

1 

2 

3 

4 

5 

6 
7 


PROBABILITY LEVEL 
0.10 
0.20 
u. 30 
0.40 
0.50 
0.60 
0.70 




JJ 

8 

9 

10 

11 

12 

13 


PROBABILITY LEVEL 
0. 80 
0.90 
0.95 
0.98 
0.99 
0.999 


CALCULATE T STATISTICS 


220 WRITE (LIST, 230) 

230 FORMAT ( 1H0 , 23HCALCUL A TED T STATISTICS /75H THE T STATISTICS CAN BE 
1 USED TO TEST THE NET REGRESSION COEFFICIENTS B(I). ) 

DO 260 J=l, NOTERM 
DO 240 K= 1 , NODEP 
TT(J,K)=ABS(B(J,K)/DEVB(J»K) ) 

240 CONTINUE 

WRITE (LIST, 250) ( TT< J » K ) ,K = 1, NODEP) 

250 FORMAT ( 1H 9G14.6) 

260 CONTINUE 


C 



C SEARCH THE TABLE OF TABULATED DEGREES OF FREEDOM 


91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 
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c 153 

MAKENJ=. FALSE. 154 

IF ( NDEG-30 >290,290,300 155 

290 I I=NDEG 156 

GO TO 4-00 157 

300 IFINDEG-4 0)310, 320, 330 158 

310 FINV=l. 0/40.0 159 

FM1INV=1. 0/30.0 160 

MAKENJ=.TRUE. 161 

320 11=31 162 

GO TO 400 163 

330 I F ( NOE G-60 1340,350,360 164 

340 FINV=1. 0/60.0 165 

FMlINV=1.0/40.0 166 

MAKEMJ=.TRUE. 167 

350 11=32 168 

GO TO 400 169 

360 IF(N DEG- 120)370, 380, 390 170 

370 FINV=1 .0/120.0 171 

FM1 I Nt/ = 1 . 0/60.0 172 

MAKENJ=.TRUE. 173 

380 11=33 174 

GO TO 400 175 

390 11=34 176 

FINV=0.0 177 

FM1IN/=1. 0/120.0 178 

MAKENU=.TRUE. 179 

180 

181 

400 WRITE(LIST,410) 182 

410 FORMAT ( 104H UNDER NULL HYPOTHESIS THE INTERVAL I-T,T) WHERE T IS G 183 

XIVEN ABOVE, HAS APPROX PROBABILITY LISTED BELOW. /42H MINUS S 184 

XIGN INDICATES PROB EXCEEDS .999. ) 185 

IF(.NDT.MAKENU) GO TO 430 186 

FNOEG= NDEG 187 

DO 420 J J=1 , 1 3 188 

T(35,JJ)=T( I I ,JJ ) + ( ( 1.0/FNDEG - FI NV ) / ( FM1 IN V-F I NV ) ) * ( T ( 1 1 -1 , J J ) 189 

1 -Till , JJ) ) 190 

420 CONTINUE 191 

11=35 192 

430 DO 550 J= 1 , NUT ERM 193 

DO 540 K= 1 , NODtP 194 

DO 440 J J= 1 , 1 3 195 

JJ=JJ 196 

I F ( T 1 1 I, JJ)-TT( J,K> ) 440,450, 460 197 

440 CONTINUE 198 

PROB( J,K)=-0.999 199 

GO TO 540 200 

450 PROBI J,K)=PLEVEL( JJ) 201 

GO TO 540 202 

460 IFIJJ.LE.9) GO TO 470 203 

JJl=JJ-2 204 

J J2= JJ -1 205 

J J3 = J J 206 

GO TO 490 207 

470 IFIJJ.LE.4) GO TO 480 208 

JJ1=JJ-1 209 

J J2 = JJ 210 

JJ3=JJ+1 211 

GO TO 490 212 

480 JJ1=JJ 213 

JJ2=JJ+1 214 

J J3= J J+2 215 
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216 

PERFORM A THREE-POINT LAGRANGE INTERPOLATION 217 

218 

490 X=AL03(TT( J,K> ) 219 

X1=AL3S(T ( I I , J J 1 ) ) 220 

X2=AL3S(T I I I , JJ2 > ) 221 

X3=AL3G(T( I I » JJ3 ) ) 222 

IF(TT( J,K).LE.1.0) GO TO 500 223 

Y1 = AL0G(1.0-PLEVEL(JJ1) ) 224 

Y2=AL3GI1.0-PL£V£L(JJ2) ) 225 

Y3=AL3G(1.0-PLEVEL( JJ3) ) 226 

GO TO 510 227 

500 Y1=AL3G(PLEVEL( JJ1) ) 228 

Y2 = AL3G(PLEVEL( JJ2) ) 229 

Y3 = AL3GIPLEVEL( JJ3) ) 230 

510 PR0B(J,K)= l (X-X2)*(X-X3)*Y1)/((X1-X2)*(X1-X3) > + ( ( X-Xl ) * ( X-X3 ) 231 

1 *Y2 ) / ( ( X2-X1 ) * ( X2-X3 ) ) + ( ( X-Xl ) * ( X-X 2 ) * Y3 ) / ( ( X3-X1 ) * ( X3-X2 ) ) 232 

IF (TT< J,K)-i.O) 520,520,530 233 

520 PROB(J»K)=EXP(PROB(J*K) ) 234 

GO TO 540 235 

530 PROB(J,K)=1.0-EXP(PROB(J»K) ) 236 

540 CONTINUE 237 

*********************************************************************** 238 

WRITE THE PROBABILITIES (1.0-ALPHA) 239 

WRITE! LIST, 550) IDOUT ( J ) , ( PROB I J, K ) , K=1 , NODEP ) 240 

550 FORMAT ( 1 H 1 3 , 9 ( 8X , F6 . 3 > ) 241 

560 CONTINUE 242 

243 

*********************************************************************** 244 

LIST THE DESIRED VALUE OF PROBABILITY ( P WANT ) 245 

246 

570 PERCEN=PWANT*100.0 247 

WR ITE(LIST,580) PERCEN 248 

580 FORMAT ( 1H0,36HTHE DESIRED VALUE OF PROBABILITY IS ,F5.1, 8H PERCEN 249 

IT ) 250 

251 

DELcTE THE TERM WITH THE LOWEST COMPUTED PROBABILITY IF THAT 252 

PROBABILITY IS LESS THAN THAT DESIRED ( P WANT ) 253 

254 

IF(. NOT. DELETE) GO TO 660 255 

IFINKEEP.EQ. NOTERM) GO TO 660 256 

I OU TO 257 

590 AMI N=PWANT 258 

JLO=MAXO( l.NTKEEP) 259 

DO 623 J=JLO, NOTERM 260 

IF(A8S(PR0B(J,1) )-PW ANT) 600, 620,620 261 

600 IF( ABS(PROB( J,l) ) -AM I N ) 6 10 , 6 20 , 620 262 

610 AMIN=4BS( PROB ( J , 1 ) ) 263 

I OU T = J 264 

620 CONTINUE 265 

IF(IOJT) 66 j , 66G , 630 266 

630 WRITE! LIST ,t>50 ) IOOUT(IOUT) 267 

650 FORMAT ( 1H 1UX.11HTHE TERM XI,I2,18H) IS BEING DELETED ) 268 

GO TO 670 269 

C ALL VARIABLES REMAINING HAVE BEEN CONCLUDED SIGNIFICANT 270 

660 RETURN1 271 

670 RETURN 272 

END 273 
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tIBFTC RSTATX 


ON 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

SUBROJTINE RSTATS ( IORDER ) 

C** ***************************** *********** 

C0MM0N/BIG/B(60, 9),SUMXY(60, 9 ), SUMXX ( 1830 ) , SUMXXI (1830) 

X ,DUMMY(1) 

OOUBLE PRECISION B, S UMXY , SUM XX , SUMXX I 

C0MM0N/MED/B0(9) , SUMX ( 69 ) , SUMX2 ( 69 ) ,SUMY2( 9) ,ZEAN(69) , 

X CONI 99) ,ERRMS(9),IDENT(13) , I OQUT ( 60 ) , NCON ( 200 ) ,NTERM(69) 
X NTRAN(IOO) , NXCOO ( 100 ) » POOL ED ( 9 ) , REP VAR ( 9 ) , RE SMS ( 9 ) , X ( 99 ) 

DOUBLE PRECISION BO, SUMX, SU MX2 , SUM Y2 , Z E AN 
COMMON /FRMTS/ FMT( 1 3 ) , FM TTR I ( 14 ) 


SUBROUTINE RSTATS 
PURPOSE 

1) COMPUTE AND PRINT THE ANALYSIS OF VARIANCE TABLES 
REGRESSION AND LACK-OF-FIT IF APPROPRIATE. 

2) COMPUTE AND PRINT R-SQUARED AND STANDARD ERROR OF 
ESTIMATE 

3) COMPUTE AND PRINT SUMS OF SQUARES DUE TO EACH VARIABLE 
IF IT WERE LAST TO ENTER REGRESSION 

4) COMPUT E AND PRINT THE STANDARD DEVIATIONS OF EACH 
REGRESSION COEFFICIENT. 


X ERRFXD, ECONMY , IOUT, ICOL 
LOGICAL ECONMY 
LOGICAL SATRTD 

LOGICAL BYPASS, BZERO, 

XIFSSR, I FTT , IFWT , 

XSTORYC, STORYX, STORYI, 

DOUBLE PRECISION RWT , TOT WT , W E I GHT 


DELETE, 

REPS, 


IFCHI, 
PREDCT , 


FIRST , ERRFXD 


DIMENSION SSQREG l 9 ) , SSQRES ( 9 ) , REGMS ( 9 ) , 

X XLOF (9 ) , XL OFMS ( 9 ) , FRATI0I9), RSQD(9), R(9), 

X SSQLSTI9), DEVB ( 60, 9 ) 

EQUIVALENCE ( DUMMY ( 1 0 ), S SQR E S ), ( DUMMY! 19 ), REGMS ) , (DUMMY! 37) , XLOF ) 

X , (DUMMY (46), XLOFMS) , ( DUMMY ( 55 ) , FR AT 10 ) , (DUMMY (64) , RSQD ) , 

X (DUMMY (73) ,R) , ( DUMM Y ( 82 ) , SSQL ST ) , ( DUMM Y ( 9 1 ) ,DE VB ) , 

X (DUMMYI700) .SSQRES) 

DOUBLE PRECISION SSQREG 

• ••**•»***»•••»*••*••*•*••»•• ** #* **»*********»************»**■*»****.*** 

COMPUTE DEGREES OF FREEDOM AND RECIPROCALS 

NREG= IORDER 

NTOT = I F I X ( TOTWT ) -1 

IFl .NOT. BZERO) NTOT = NTOT + 1 

NRE S= NTOT-NREG 

NLUF= NRES - NPDEG 

RNREG= 1 . O/FLOAT ( NREG ) 

IF< NRES.EQ.G ) GO TO 980 
RNRES= 1.0/FLUAT (NRES ) 

SATRT3=. FALSE. 

IFCjLDF.FQ.O) GO TO 90 


10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
23 



COMMON/SMALL/ 

BYPASS, 

BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

24 

X 

IFTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

25 

X 

ISTRAT, 

JCOL , 

KONNO, 

LENGTH, 

LIST, 

26 

X 

NERROR, 

NODEP, 

NOOB, 


NOTERM, 

27 

X 

NOVAR, 

NPDEG, 

NRES, 

NTRANS , 

NWHERE , 

28 

X 

P, 

PREDCT, 

REPS, 

RWT, 


29 

X 

STORYI , 

STORYC 

, STORYX, 

TOTWT, 

WEIGHT, 

30 


31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 
5 7 

58 

59 

60 


60 
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RNLUF=1.0/FLOAT(NLOF) 61 

GO TO 100 62 

90 SATRTD=.TRUE. 63 

100 CONTINUE 64 

NXTFRM=I0RD£R 65 

RN00B= RWT 66 

67 

68 

COMP J TE RESIDUAL SUM OF SQUARES, RESIDUAL VARIANCE, VARIANCE 69 

FROM REPLICATIONS IF APPROPRIATE, AND THE F-RATIO OF MEAN SQUARE 70 

LACK-3F-F I T AND MEAN SQUARE RESIDUALS. 71 

DO 210 J= 1 , NODEP 72 

SSQRE3 ( J ) =0. 0 73 

DO 200 1=1, NXTERM 74 

SSQRESl J)=SSQREGl J) + B ( I, J ) *SUMXY ( I , J ) 75 

200 CONTINUE 76 

SSQRES (J)=SUMY2(J)-SSQREG( J ) 77 

REGMS ( J ) = SSQREG(J)* RNREG 78 

RESMSl J) = SSQRESI J)*RNRES 79 

RSQD( JJ=SSQREG(J )/SUMY2( J) 80 

R(J)=SQRT (RSQD(J) ) 81 

IF( ( .NOT. RFPS ) .OR.SATRTD ) GO TO 210 82 

XLOF( J)=SSQRES( J)-POOLED(J) 83 

XLOFMS ( J ) = XLOFI J)*RNLOF 84 

FRATI31 J)=XLOFMSl JI/REPVARl J ) 85 

210 CONTINJE 86 

87 

•••if******************************************************************* 88 

DETERMINE WHICH ESTIMATE OF SIGMA SQUARED SHOULD BE USED IN 89 

HYPOTHESIS TESTS. PUT THE PROPER ONE IN ERRMS AND SET ERRFXD 90 

TO TRUE IF THE PRESENT VALUE IS TO BE USED FOR ALL FOLLOWING 91 

TESTS AND T-STATIST ICS. 92 

I OUT= I STR AT 93 

IF ( ERRFXD) GO TO 250 94 

IFl ISTRAT.NE.3I GO TO 214 95 

211 DO 213 J=l, NODEP 96 

213 ERRMS ( J ) = RESMS(J) 97 

NERROR = NRES 98 

I OUT = 3 99 

GO TO 250 100 

214 IFl ISTRAT.NE.l) GO TO 218 101 

IFl .NOT. REPS) GO TO 211 102 

00 215 J=l, NODEP 103 

215 ERRMS l J ) = REPVARIJ) 104 

NERROR = NPOtG 105 

ERRFXD= .TRUE. 106 

1 OUT = 1 107 

GO TO 250 108 

218 IFl FIRST. AND. I I ORDER . EQ. NOTERM ) ) GO TO 220 109 

GO TO 211 110 

220 ERRFX0= .TRUE. Ill 

DO 222 J=l, NODEP 112 

222 ERRMS l J ) = RE SMS l J ) 113 

NERROR = NRES 114 

I STRAT = 2 115 

I 0UT=2 116 

117 

118 

WRITE ANOVA TABLES 119 

250 DO SOD J=l, NODEP 120 

IFIERRMSI Jl.EQ.O.O) ERRMSl J ) =1. 0E-30 121 

WRITEl LIST, 1001 ) J 122 
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WRITEILIST, 1002) 123 

WRI TE( LIST, i003) SSOREG(J), NREG, REGMSIU) 124 

WRITE! LIST, 1004) SSQRES(J), NRES » RESMSIJ) 125 

WRITE! LIST, 1005) 126 

WRITE! LIST, 1006) SUM Y2 < J ) , NTOT 127 

WRI TE1 LIST , 1007) 128 

WRITEILIST, 1500) RSQO(J), RIJ) 129 

ST0=S3RTIRESMS! J) ) 130 

WRITEILIST, 1600) STD 131 

WRITE! LIST, 1700) IOUT , ERRMS ! J ) , NERRDR 132 

F=REGMSl J l/ERRMS! J) 133 

WRITE! L I ST, 1 750 IF, NR EG, N ERROR 134 

IF! I.NOT.REPSI.OR.SATRTD) GO TO 500 135 

WRITE! LIST, 2001) 136 

WRITE! LIST, 1002) 137 

WRITE! LIST, 2005) XLOF!J), NLOF , XLOFMSIJ) 138 

WRITE! LIST, 2006) POOLED!J), NPDEG, REPVARIJ) 139 

WRITEILIST, 1004) SS QRES ! J ) , NRES , RE SMS ( J ) 140 

WRITE! LIST, 1005) 141 

WRITE! LIST, 2008) FRATIO(J) 142 

WRITEILIST, 1007) 143 

500 CONTINUE 144 

145 

a-***#****#*************#**-********#************************************ 146 

COMPUTE CONTRIBUTION OF EACH INDEPENDENT VARIABLE TO REG SUM 147 

OF SQJARES AS IF IT WERE LAST TO ENTER 148 

WRITE! LIST, 370) 149 

I R = 0 150 

DO 8635 K=1 , NXTERM 151 

IR= IR+K 152 

DO 3632 J = 1 , NOOEP 153 

8632 SSQLSI ( J ) = 8 I K , J ) **2 / SUMXX I ! IR ) 154 

WRI TE( LIST, 380) IDOUT l K ) , < SSQLST { J ) , J = 1 , NODEP ) 155 

8635 CONTINUE 156 

157 

*********************************************************************** 158 

COMPUTE STANDARD DEVIATION OF REGRESSION COEFFICIENTS 159 

WRITE! LIST, 375) 160 

IF! , NOT. 8ZER0 ) GO TO 959 161 

DO 913 J=l, NXTERM 162 

R ( J ) =0 « 0 163 

DO 913 1=1, NXTERM 164 

CALL LOC! I , J , I R ) 165 

R ( J ) = 1 1 J ) + ZE AN ( I ) *S'J MXX I { I R ) 166 

910 CONTINUE 167 

XXT = 0,3 168 

DO 923 J=l, NXTERM 169 

920 XXT=XXT+ZEANI J)*RI J) 170 

DO 933 K=l, NODEP 171 

930 DEVB(l,K)=SQRTl ERRMS l K ) * { RNOOB+XXT ) ) 172 

K=0 173 

WRITE! LIST, 380) K ,! D EVB < 1 , J ), J = 1 , NODEP ) 174 

959 IR=0 175 

DO 973 J=l, NXTERM 176 

I R= IR+J 177 

DO 963 K=l, NODEP 178 

DEV B ( J , K ) =SQRT( ERRMS(K) *SUMXXI l IR ) ) 179 

960 CONTINUE 180 

WRITEILIST, 380) I DOUT ( U ) , I D EVB ( J , KR ) , KR = 1 , NODEP ) 181 

970 CONTINUE 182 

RETURN 183 

C 184 
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*********** ************************************************** ********** 185 

FORMATS 186 

1001 FORMAT (42H4AN0VA OF REGRESSION ON DEPENDENT VARIABLE 15) 187 

1002 FORMAT ( 1H 79(1H*)/79H SOURCE SUMS OF SQUARES DEG 188 

XREES 3F FREEDOM MEAN SQUARES /1H 79(1H-)) 189 

1003 FORMAT C 17H REGRESSION G20.8, 5X , 1 10 , 5 X , G20. 8 ) 190 

1004 FORMAT ( 17H RESIDUAL G20.8, 5X, 1 10 , 5X , G2 0. 8 ) 191 

1005 FORMAT ( 1H 79(lH->) 192 

1006 FORMAT ( 17H TOTAL G20.8» 5X.U0) 193 

1007 FORMAT ( 1H 79(1H*)) 194 

2001 FORMAT (1X/1X/22H ANOVA OF LACK OF FIT ) 195 

2005 FORMAT ( 17H LACK OF FIT G20.8, 5X , 1 10 , 5X, G2 0. 8 ) 196 

2006 FORMAT ( 17H REPLICATION G20.8, 5X , I 10 , 5X, G2 0. 8 ) 197 

2008 FORMAT ( 28H F = MS ( LOFJ /MSI REPS ) = F10.3 ) 198 

370 FORMAT (74H1 SUMS OF SQUARES DUE TO EACH VARIABLE IF IT WERE LAST T 199 

XO ENTER REGRESSION ) 200 

375 FORMAT ( 1 15H2 STANDARD DEVIATION OF REGRESSION COEFFICIENTS (DERIVE 201 

XD FROM DIAGONAL ELEMENTS OF (X TRANSPOSE X) INVERSE MATRIX )) 202 

380 FORMAT ( 1H I3.9G14.6) 203 

1500 FORMAT ( 40H R SQUARED = SSQ(REG) / SSQ(TOT) = F8.6, 204 

X 5X, 4-HR = F7.6) 205 

1600 FORMAT ( 34H STANDARD ERROR OF ESTIMATE G14.6) 206 

1700 FORMAT ( 24H USING POOLING STRATEGY I2.25H THE ERROR MEAN SQUARE * 207 

X G14.7. 26H WITH DEGREES OF FREEDOM = 16) 208 

1750 FORMAT (5X,19HF=MS(REG)/MS(ERR)= F 6. 2 , 5X , 1 3HC0MP ARE TO F(I2,1H,I3,1 209 

XH)) 210 

980 WRITE(LIST,981) 211 

981 FORMAT ( 41H ZERO RESIDUAL DEGREES OF FREEDOM. STOP. ) 212 

STOP 213 

END 214 


ilBFTC RECTX 


SUBROUTINE RECT ( I ROW , J JC OL , I MAX , JM AX , A , B . FM T , I I > 1 

DIMENSION A( IMAX, JMAX) , FMT( 14) ,XOUT( 8) 2 

DOUBLE PRECISION B.DXOUT 3 

DIMENSION B( IMAX, JMAX ) .DXOUT ( 8 ) 4 

COMMON/SMALL/DUM( 15) .LIST 5 

DATA J 8/ 8/ 6 

LOGICAL OUT 7 

OUT =. FALSE. 8 

JT I MES =0 9 

JCOL= J JCOL 10 

5 JNXT = J COL— J8 U 

IF(JNXT) 10,20,30 12 

10 JP= JCOL 13 

GO TO 40 14 

20 JP=J8 15 

GO TO 40 16 

30 JCOL=JNXT 17 

JP=J8 18 

GO TO 50 19 

40 OUT = . TRUE. 20 

50 DO 100 1 = 1,1 ROW 21 

GO TO (55,75) , I I 22 

55 CONTINJE 23 

DO 60 J=1,JP 24 

J J = J TI MES +J ?5 
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60 XQUHJ) = A ( I , J U ) 

WRI TE( LIST.FMT) I , ( XOUT ( K ) ,K=1, JP) 

go rn 100 

75 DO 50 U=1,UP 
J J = J T I MES + J 
80 OX OUT ( J ) = B ( I , JU ) 

WRI TE( LIST ,FMT ) I, (DXOUTtJ ), U = l, JP) 
100 CONTINUE 

I F ( OUT ) RETURN 
WR I T E ( L 1ST, i 10 ) 

110 FOR M AT ( 1H /1H ) 

UTIMES=JTIMES +JP 

GO TO 5 

END 


26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 
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C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE PRtDIC 
PURPOSE 

DREAD INPUT LEVELS OF INDEPENDENT VARIABLES AND COMPUTE 

A PREDICTED RESPONSE FROM THE ESTIMATED REGRESSION EOUATIUN. 
2) COMPUTE VARIANCE AND STADARD DEVIATION OF THE PREDICTED 
MEAN VALUE AND A SINGLE FURTHER OBSERVATION. 

SUBROUTINES NEEDED 
TRANS 
LOC 

REMARKS 

VALUES FOR DEPENDENT VARIABLES ARE NOT NECCESARY FOR THE 
PREDICTING OF VALUES. HOWEVER, A DUMMY VALUE MAY NEED TO 
BE SUPPLIED IF A ZERO (BLANK) INPUT VALUE WILL CAUSE AN 
IMPOSSIBLE OPERATION TO BE ATTEMPTED DURING THE 
TRANSFORMATIONS. 

,*************************************»***************»**»**»»»»»»»* 

SUBROUTINE PREDIC 
COMMON /FRMTS/ FMT ( 1 3 ) , FMTTRI ( 14 ) 

COMMOV/BIG/B(6G,9 ) , S U MX Y ( 60 , 9 ) , SUMXX ( 18 30) , SUMXXI ( 1830) 

DOUBLE PRECISION B , SUMXY , SUMXX , SUMXX I 

COMMUN /MED/801 9 ) , SUM X ( 69 ) , SUMX 2 I 69 ) .SUMY2I9) ,ZEAN(69) , 

X CONI 99) , ERRMS ( 9 ) , I DENT (13) , IDOUTI &0) , NCUN(200) , NTERM169) , 

X NTRANt 100) ,NXCOD( 100) , POOL E D ( 9 ) , R C P V AR ( 9 ) , R E SMS ( 9 ) , X ( 99 ) 


X ERKFXD, ECONMY, 
LOGICAL ECONMY 
LOGICAL BYPASS, 

XIFSSR, IFTT, 

XSTORYC , STORYX, 


I OUT , ICOL 

<ZERO, 
IFWT, 
STORYI , 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
2 7 
28 



DOUBLE PRECISION 

BO, SUMX, SU MX2,SUMY2, 

ZE AN 


29 


COMMON/SMALL/ 

BYPASS ,BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

30 

X 

IFTT, 

IFWT, 

I NPUT, 

INPUTS, 

INTER, 

31 

X 

ISTRAT, 

UCOL , 

KONNO , 

LENGTH , 

LIST, 

32 

X 

NERROR, 

NODCP, 

NOOB, 


NO TERM, 

33 

X 

NOVAR, 

NPDEG, 

NRES, 

NTRANS , 

4 si HERE , 

34 

X 

P, 

PREDCT, 

REPS, 

RWT, 


35 

X 

STORYI , 

S TORYC, 

STORYX, 

r o f w t , 

W rr I GHT , 

36 


DELETE, IECHI, 

REPS, PREUCT, 

FIRST , ERR E X 0 


3 7 
38 
3 9 

40 

41 
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l 


DOUBLE PRECISION WEIGHT, RWT, TOTWT 42 

DIMENSION YC ALC ( 9 ) , V(60), VARM(9), SEEM(9), 43 

X V ARP ( 9 ) , SEEPI9) 44 

DOUBLE PRECISION XXT , V 45 

EQUIVALENCE I YC ALC I 1 ) , SUMXX (1) ) , ( V (1) , SUMXX { 1 50 ) ) , 46 

X ( V ARM ( 1) ,SUMXX(71 ) ) , ( SEEM ( i ), SUMXX ( 80 )) , ( VARP ( 1 ), SUMXX ( 89 ) ) 47 

X , ( SEEP (1 ) ,SUMXX(98) ) 48 

EQUIVALENCE (RNOOB.RWT) 49 

50 

52 

IF{ P10TERM. EQ.O) RETURN 53 

WRI TE C LIST, 3 ) 54 

READ ( I NPUT5 , 5 ) NPRED 55 

56 

DO 500 KK=1, NPRED 57 

58 

105 READ! I NPUT5, FMT ) (XII) , 1 = 1, ICOL ) 59 

WRITE! LIST, 110) (XII) , 1 = 1, ICOL) 60 

125 CALL TRANS 61 

DO 130 K= 1 , J COL 62 

I=NTERM( K ) 63 

X ( K ) = CON (I) 64 

130 CONTINUE 65 

WRI I E ( LIST, 135) t X ( I) , 1 = 1, NOTERM) 6 6 

67 

COMPUTE PREDICTED RESPONSE 68 

140 DO 150 K= 1 , NODEP 69 

YC ALC ( K ) = 80 ( K ) 70 

I F ( . NOT. BZERO ) YCALC ( K ) = .j»0 71 

DO i5 0 J=l, NOTERM 72 

YCALC ( K ) = YCALC ( K ) + B(J,K)*X(J) 73 

150 CONTINUE 7 4 

75 

COMPUTE VARIANCE AND STANDARD DEVIATION OF REGRESSION LINE 76 

AND VARIANCE AND STANDARD DEVIATION OF PREDICTED VALUE 77 

AT THE POINT XO 78 

79 

DO 250 K=l, NOTERM 80 

V(K)=O.ODO 81 

DO 250 J=l, NOTERM 82 

CALL L OC ( J , K » I R ) 83 

V ( K ) =V (K) + (XI J )-ZEAN{ J ) )*SUMXXI ( IR) 84 

250 CONTINUE 85 

XXT=0. 0D0 86 

DO 275 K=l, NOTERM 87 

XXT = XXT + (X(K)-ZEAN(K)) *V(K) 88 

275 CONTINUE 89 

X RNUOB = RNOOB 90 

IF(. NOT. BZERO) XRN008=0.0 91 

DO 300 K=l, NODEP 92 

VARM ( < ) = ERRMS ( K ) *(XRNOOB + XXT) 93 

SEEM ( X ) =SQRT ( VARM ( K ) ) 94 

VARP ( < ) = ERRMS(K)+VARM(K ) 95 

SEEP(<)=SQRT(VARP(K) ) 96 

300 CONTINUE 97 

WRI TE( LI ST, 3 10) ( YCALC (K ) ,K = 1 , NODEP ) 9 8 

WRITE! LIST, 320) (VARM ( K ), K = 1 , NODEP ) 99 

WRI TE! LIST, 320 ) I SEEM ( K ), K = 1 , NODEP ) 100 

WRI TE(LIST, 320) (VARP ( K ), K = 1 , NODEP ) 101 

WRITE! LIST, ^20) (SEEP ( K ) ,K = 1 , NODEP ) 102 

C 103 
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c 


104 

500 CONTINUE 105 

RETURN 106 


3 FORMAT ( 54H1F0R EACH SET OF INDEP VARIABLES THERE IS COMPUTED... / 


X 2 OH PREDICTED RESPONSE / 
X 29H VARIANCE OF REGRESSION LINE / 
X 34H STANDARD DEVIATION OF REGRESSION / 
X 29H VARIANCE OF PREDICTED VALUE / 
X 39H STANDARD DEVIATION OF PREDICTED VALUE ) 

5 FORMAT ( 14 ) 


110 FORMAT (39HKINPUT DATA FOR THIS PREDICTED RESPONSE /(1H 9G14.6)) 
135 FORMAT ( 56HK INDEPENDENT TERMS ACCORDING TO FINAL REGRESSION MODEL 
X /(H 9G14.6) ) 

310 FORMAT! 55HKPREDICTED RESPONSE FOR ABOVE INDEP VARIABLES 
X /<1H 9G14.6 ) ) 

320 FORMAT ( 1H 9G14.6) 

END 


107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 


* I BFTC SUMUPX 


SUBROUTINE SUMUPS 
PURPOSE 

1) CALCULATE (X TRANSPOSE X) AND (X TRANSPOSE Y) MATRICES ONE 
OBSERVATION AT A TIME. 

2) COMPUTE TOTAL OF THE WEIGHTS 

** BOTH CALCULATIONS ARE IN DOUBLE PRECISION 


SUBROUTINES 

LOC 


NEEDED 


SUBROUTINE SUMUP 

COMMO\I/BIG/S(60,9) » S UMX Y (60*9) * SUM XX ( 1830) , SUMXX I ( 1 830 ) 
DOUBLE PRECISION B , S UMXY , SUM XX, SUMXX I 

COMMON/MED/BO! 9 ) ,SUMX ( 69 ) , SUMX2 ( 69 ) , SUMY2 ( 9 ) , ZEAN ( 69 ) , 

X CON! 99) , ERRMS ( 9 ) , I DENT (13) , I DOUT { 60 ) , NCON ( 200 ) , N T E R M ( 6 9 ) , 
X N TRAN! 100 ) ,NXCOD( 100 ) , POOL ED ( 9 ) , REP VAR ( 9 ) , RtSMS ( 9 ) , X ( 99 ) 
DOUBLE PRECISION BO,SUMX,SU MX2 , SUM Y2 , Z E AN 


LOGICAL 
LOGICAL 
XIFSSR, 
XSTORYC , 


ECONMY 

BYPASS, 
I FTT , 
STORYX, 


BZERO, 
IFWT , 
STORYI , 


DELETE, 

REPS, 


IFCHI , 

PREDCT 


FIRST , ERRFXD 


DOUBLE PRECISION RWT , TOT WT , WE I GHT 
DOUBLE PRECISION DUB1,DUB2 


DO I ID 1=1, JCOL 

SUMX ( I ) =SUMX ! I )+X ( I ) * WE I GHT 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


COMMON/SMALL/ 

BYPASS, 

BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

21 

X 

IFTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

22 

X 

ISTRAT, 

UCOL, 

KONNO, 

LENGTH, 

LIST, 

23 

X 

NERROR, 

NODEP, 

NOOB, 


NOTERM, 

24 

X 

NOVAR, 

NPDEG, 

NRES, 

NTRANS , 

NWHERE , 

25 

X 

P, 

PREDCT, 

REPS, 

RWT, 


26 

X 

STORYI , 

STORYC 

, STORYX, 

TOTWT , 

WEIGHT, 

27 

X 

ERRFXD, ECONMY, 

I OUT , 

ICOL 



28 


29 

30 

31 

32 

33 

34 

35 

36 

37 

38 
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on on 


110 CONTINUE 3SJ 

I R=0 40 

DO 103 K=l, NOTERM 41 

DUB1=XU4J 42 

00 90 J=l,NODEP 43 

KBAR=J +N0TERM 44 

DUB2=X(K8AR) 45 

SUMX Y{ X, J )=SUMXY(K, J > +DU B1 *DUB2 * WE I GHT 46 

90 CONTINUE 47 

DO 50 1=1, K 48 

C 49 

I R= l R+ 1 50 

DUB2=X ( I ) 51 

SUMXX( I R ) =SUMXX ( I R H- DUB1 *DUB2*WE I GHT 52 

50 CONTINUE 53 

100 CONTINUE 54 

00 15 J = 1 » NODEP 55 

KBAR=NOTERM + J 56 

DUB1=X ( KBAR ) 57 

SUMY2I J)=SUMY2( J)+DUB1*DUB 1* WEIGHT 58 

15 CONTINUE 59 

TOT WT= TOTWT + WE I GHT 60 

RETURN 61 

END 62 


$IBFTC BORDXX 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


SUBROUTINE BORD 
PURPOSE 

TO COMPLETE THE INVERSION OF A SYMMETRIC POSITIVE DEFINITE 
MATRIX A OF ORDER N GIVEN THAT THE UPPER LEFT SUB- 
MATRIX OF ORDER N-l HAS ALREADY BEEN INVERTED. 

SUBROUTINES NEEDED 
LOC 


REMARKS 

ONLY THE UPPER TRIANGULAR PART OF A IS STORED AS A 
VECTOR IN THE ORDER A ( 1 , 1 ) , A ( 1 , 2 ) , A I 2 , 2 ) , A ( 1 , 3 ) , . . 
SUBROUTINE BORD I I ORD ER , A ) 


ETC 


DIMENSION BETA! 60), All) 

DOUBLE PRECISION A, AL PHA , R AL PHA .BETA 


ALPHA= 0. 000 
NM1 = I ORDER-1 
IF(NMl) 100,100,200 
100 A ( 1 ) = 1.0/AU) 

GO TO 600 

200 M=NMl»(NMl+i)/2 
LEN = M + I ORDER 
C 

DO 403 1=1, NM1 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 
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BETACI )= 0.0 DO 31 

M I = M+I 32 

DO 350 J= 1 . NMi 33 

CALL LOCC I, J, II ) 34 

MJ= M«-J 35 

BETA { I ) = BETA ( I )-A( I I )*A(MJ) 36 

350 CONT I N JE 37 

ALPHA= ALPHA + A ( M I ) *BET AC I ) 38 

400 CONTINUE 39 

40 

41 

ALPHA = ALPHA + ACLEN) 42 

RALPHA=1.0 DO /ALPHA 43 

ACLEN) = RALPHA 44 

45 

DO 500 1=1, NMI 46 

DO 500 J=1,I 47 

CALL LOCC I,J,II) 48 

A(II)= A ( I I ) + BETAC I ) *BETA( J)*RALPHA 49 

500 CONTINUE 50 

51 

DO 550 J= 1 , NM 1 52 

MJ= MO 53 

A ( M J ) = BETAC J)*RALPHA 54 

550 CONTINUE 55 

56 

57 

600 CONT I N JE 58 

RETURN 59 

END 60 


IBFTC MFIXXX 

SUBROUTINE MFIX 1 

2 

3 

COMMON /FRMTS/ FMT ( 1 3 I , FMTT R I ( 14 ) 4 

COMMON/BIG/B (60,9) ,SUMXYC60, 9) , SUMXXC 1830) , SUMXXI ( 1830) 5 

DOUBLE PRECISION B, S U MX Y , SUMXX , SUMXX I 6 

COMMON /MED/ BOC 9) , SUMX ( 69 ) , SUMX2 C 69 ) , SUMY2 C 9 ) , ZEAN ( 69 ) , 7 

X CONC99) ,ERRMS(9),IDENT(13) , I DOUT ( 60 ) , NCON ( 200 ) ,NTERM(69) , 8 

X NTRANC100) , NXCOD ( 1 OC ) , POOL ED ( 9 ) , RE PV AR ( 9 ) , RE SMS l 9 ) , X ( 99 ) 9 

DOUBLE PRECISION BO, SUMX, SU MX2 , SUM Y2 , Z E AN 10 



COMMON/SMALL/ 

BYPASS, 

BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

11 

X 

IFTT, 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

12 

X 

ISTRAT, 

JCOL, 

KONNO, 

LENGTH, 

LIST, 

13 

X 

NERROR, 

NODEP, 

NOOB , 


NOTERM, 

14 

X 

NOVAR, 

NPDEG, 

NRES , 

NTRANS, 

NWHERE, 

15 

X 

P, 

PREDCT , 

REPS, 

RWT, 


16 

X 

STORYI , 

STORYC 

, STORYX, 

TOT WT , 

WEIGHT, 

17 


X ERRFXD, ECONMY , IOUT, ICOL 
LOGICAL ECONMY 

DOUBLE PRECISION RWT , TOT WT , W E I GHT 
LOGICAL BYPASS, HZERO, 

XIFSSR, IFTT , IFWT , 

XSTORYO, STORYX, STORYI, 


DELETE, 

REPS, 


IFCHI , 

PREOCT, 


FIRST , ERRFXD 




18 

19 

20 
21 
22 

23 

24 
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C 25 

IF ( EC3NMY ) GO TO 500 26 

IF(GZERO) GO TO 500 27 

WRITE! LIST, 530) 28 

WRITE (LIST, 540) (SUMX ( I ) » I = 1 » JCOL ) 29 

WRI TE( LIST, 560 ) 30 

CALL TRIANG(X,SUMXX»NOTERM,8,FMTTRI,2) 31 

WRI TE( LIST, 565) 32 

CALL RECT(N0TERM,N0DEP,60,9,X,SUMXY,FMTTRI,2) 33 

*********************************************************************** 34 

COMPUTE AND PRINT MEANS. COMPUTE AND PRINT THE ( X TRANSPOSE X) 35 

MATRIX IN TERMS OF DEVIATIONS FROM MEAN. THE DEVIATIONS FORM 36 

OF (X T X) IS THE VARIANCE-COVARIANCE MATRIX OF THE 37 

INDEPENDENT VARIABLES. 38 

500 CONTINUE 39 

RWT=1. ODO/TUTWT 40 

DO 573 1=1, JCOL 41 

570 ZEAN ( I ) =SUMX ( I ) *RWT 42 

WRITE! LIST, 580) 43 

WRITE! LIST, 540) ( ZEAN (I) , I = 1 , JCOL ) 44 

IR = 3 45 

DO 603 J=l, NOTERM 46 

I R= I R + J 47 

IF! .NOT. BZERO) GO TO 601 48 

SUMX2! J)=SUMXX! IR)-SUMX( J)**2 *RWT 49 

GO TO 600 50 

601 SUMX2! J)=SUMXX( IR) 51 

600 CONTINUE 52 

602 CONTINUE 53 

I R= i 54 

DO 623 J=l, NOTERM 55 

DO 618 K= 1 , NODEP 56 

IF ( . NOT • BZERO) GO TO 618 57 

KBAR=NOTERM+K 58 

SUMXY! J,K)=SUMXY< J,K ) -SUMX ( J )*SUMX(KBAR ) *RWT 59 

618 CONTINUE 60 

619 DO 623 K= 1 , J 61 

IF! .NOT. BZERO) GO TO 6191 62 

SUMXX! I R ) = S UMXX ( I R ) - SUMX ( K ) * SUMX ! J ) *RWT 6 3 

6191 SUMXXI ( IR )=SUMXX( IR) /DSQRT ( SUMX2! J)*SUMX2(K ) ) 64 

6194 IR=IRH 65 

620 CONTINUE 66 

IF!. N3T. BZERO) GO TO 6220 67 

DO 6210 J=l, NODEP 68 

K=NOTE RM+ J 69 

SUMY2 ( J)=SUMY2( J ) -S UMX l K ) ** 2*R WT 70 

6210 CONTINUE 71 

6220 CONTINUE 72 

C*** *************** ******** ******* ************************************** 73 

IF! ECONMY) GO TO 622 74 

IF! .N3T. BZERO) GO TO 621 75 

WRI TE( LIST, 625) 76 

CALL TRIANGJX, SUMXX, NOTERM, 8, FMTTRI, 2) 77 

WRITE(LIST,630) 78 

CALL RECTINOTERM, NODEP, 6C, 9 , X , SUMX Y , FMTTR I ,2) 79 

621 WRITE! LIST, 670) 80 

CALL TRIANGIX, SUMXXI, NOTERM, 8, FMTTRI, 2) 81 

622 CONTINUE 82 

C************************** ***** **************************************** 83 

RETURN 84 

C 85 

530 FORMAT ! 1H0 , 3 2H SUMS OF INDEP AND DEP VARIABLES ) 86 
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540 FORMAT ( 1H 8G15.7) 87 

560 FORMAT ( 21H2X TRANSPOSE X MATRIX » 88 

565 FOR MAT (2 1 H2X TRANSPOSE Y MATRIX > 89 

580 FORMAT ( 33H MEANS OF INDEP AND DEP VARIABLES ) 90 

625 FORMAT ( 53H2X TRANSPOSE X MATRIX WHERE X IS DEVIATION FROM MEAN ) 91 

630 FORMAT ( 60H2X TRANSPOSE Y MATRIX WHERE X AND Y ARE DEVIATIONS FROM 92 

XMEAN ) 93 

670 FORMAT ( 25H2CORRELATI ON COEFFICIENTS ) 94 

END 95 


SIBFTC LOCXXX 


SUBROUTINE LOCCI,J»IR) 1 

I X= I 2 

JX= J 3 

20 IF(fX-JX) 22,24,24 4 

22 I RX= IX + (JX*JX-JX)/2 5 

GO TO 36 6 

24 IRX= JX + <IX*IX - 1X1/2 7 

36 I R= I R X 8 

RETURN 9 

END 10 


$ I BFTC XTRANS 

SUBROUTINE TRANS 

COMMON/ B I G7B( 60, 9 ) ,SUMXY (60, 9} ,SUMXX( 1830) , SUMXXI (1830) 
DOUBLE PRECISION B , SUMXY , SUMXX, SUMXX I 

COMMON/MED/BO( 9) ,SUMX (69 ), SUMX2I 69) ,SUMY2( 9) , ZEAN( 69) , 

X CONI 99 ) , ERRMS t 9) , I DENT! 13) , ID0UTI60) ,NCON ( 200 ), NTF RM ( 6 9 ) , 
X NT RAN ( 100) , NXCODt 100) , POOLED) 9),REPVAR(9),RESMS(9),X(99) 
DOUBLE PRECISION BO,SUMX,SU MX2 , SUMY2 , ZE AN 


X ERRFXD, ECONMY , IOUT, ICOL 
LOGICAL ECONMY 

LOGICAL BYPASS, C.ZERO, 

XIFSSR, IFTT, IFWT , 

XSTORYC, STORYX, STORYI, 

DOUBLE PRECISION RWT , TOT WT , W E I GMT 


DEL FTE , IFCHI 

REPS, PR E DC ’ 

FIRST , ERRl-XD 


***************************»*****************»**-***********#*******»*#* 

THIS SUBROUTINE PERFORMS TRANS FORMAT I UN S IF THIS OPTION IS 
REOUESTED. 


1 

2 

3 

4 

5 

6 

7 

8 
9 



COMMON/SMALL/ 

BYPASS , 

BZERO, DELETE, 

FIRST, IFCHI, 

IFSSR, 

10 

X 

IFTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

11 

X 

ISTRAT , 

JCOL , 

KONNO, 

LENGTH, 

LIST, 

12 

X 

NERROR, 

NODEP , 

NOOB, 


NOTERM , 

13 

X 

NOVAR, 

NPDEG, 

NRES, 

N TRANS, 

NWHERE , 

14 

X 

P, 

PREDCT , 

REPS, 

RWT, 


15 

X 

STORYI , 

STORYC 

, STORYX, 

TOT WT , 

WEIGHT, 

16 


17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27 
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28 

K TRANSFORMATION SET NUMBER. 29 

NCON ( 2 *K— 1 ) CONSTANT NUMBER TO USE. 30 

NCUN ( 2 *K ) DERIVED CONSTANT. 31 

NTRANIO NUMBER OF TRANSFORMATION REQUESTED. 32 

NXCOD(K) VARIABLE NUMBER 33 

34 

80 DO 503 K= 1 * NTRANS 35 

I =NCON ( 2*K-1 ) 36 

IF ( I ) 1 0 0 » 100,110 37 

100 CON S = D . 0 38 

GO TO 120 39 

no cons=:on( 1 1 40 

120 I =NXCO D ( K ) 41 

Y = X ( I ) 42 

MTRAN = NTRAN(K) 43 

IF ( MTR AN. LE. 0 ) MTRAN = 32 44 

140 GO T0( 150, 160, 170, 180, 190, 200, 210,220,230,240,250,260,270,280,290, 45 

X 300, 3 10, 320, 330, 340, 3 50, 36 0, 370,38 0,390,400,410,420,430,440, 46 

X 442.450) , MTRAN 47 

150 CONS=Y+CONS 48 

GO TO 460 49 

160 CONS=Y*CONS 50 

GO TO 460 51 

170 CONS=CONS/Y 52 

GO TO 460 53 

180 CONS = E XP ( Y ) 54 

GO TO 460 55 

190 CONS = Y **CONS 56 

GO 10 460 57 

200 CONS=ALOG(Y) 58 

GO TO 460 59 

210 C0NS=AL0G10(Y ) 60 

GO TO 460 61 

220 CONS = S I N ( Y ) 62 

GO TO 460 63 

230 CONS=COS(Y) 64 

GO TO 460 65 

240 C0NS=SIN(3.14159265*(C0NS*Y) ) 66 

GO TO 460 67 

250 CONS = :OS( 3.14159265* (CONS*Y ) ) 68 

GO TO 460 69 

260 CON S = 1 . 0/ Y 70 

GO TO 460 71 

270 CONS = EXP( CONS/Y ) 72 

GO TO 460 73 

280 CONS=EXP(CONS/(Y*Y) ) 74 

GO TO 460 75 

290 CONS = S QRT ( Y ) 76 

GO TO 460 77 

300 C0NS=1.0/SQRT(Y) 78 

GO TO 460 79 

310 CONS=CONS**Y 80 

GO TO 460 81 

320 CONS=lO.O**Y 82 

GO TO 460 83 

330 CONS=S I NH { Y ) 84 

GO TO 460 85 

340 CONS = C OSH I Y ) 86 

GO TO 460 87 

350 CONS={ l.O-COS(Y) 1/2.0 88 

GO TO 460 89 
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360 C0NS = A T AN I V ) 90 

GO TO 460 91 

370 C0NS=ATAN2 ( Y.CONS ) 92 

GO TO 460 93 

380 CONS=Y*Y 94 

GO TO 460 95 

390 CONS=Y*Y*Y 96 

GO TO 460 97 

400 CONS=ARSIN(SQRT (Y) ) 98 

GO TO 460 99 

410 C0NS=2.0*3.14159265*Y 100 

GO TO 460 101 

420 C0NS=1. 0/(2. 0*3. 14159265*Y) 102 

GO TO 460 103 

430 CON S = E RF ( Y ) 104 

GO TO 460 105 

440 C0NS=3 AMMA ( Y ) 106 

GO TO 460 107 

442 CONS=Y /CONS 108 

GO TO 460 109 

450 CONS=Y 110 

460 I=NC0N(2*K) 111 

480 CON( I ) =CONS 112 

I F ( I -SO » 500,500,490 113 

490 X(I)=;ONS 114 

500 CONTINUE 115 

RETURN 116 

END 117 


$ I BFTC OUTPLX 


************************* •*••••••**•»•••**•****••••»•**•*•*••••*****••• 

*•••••***•*••••*••*••••*••»•• ****************************************** 

SUBROUTINE OUTPLT(PNCH) 

COMMON /B I G/B(60,9),SUMXY(60,9) , SUMXXt 1830) ,SUMXXI (1830) 

DOUBLE PRECISION B, SUMXY ,SUMXX. SUMXX I 

COMMON/MED/BOl 9 ), SUMX( 69 ),SUMX2(69) , SUMY2 ( 9 ) , ZE AN ( 69 ) , 

X CONI 99) ,ERRMS(9) , I DENT (13) ,ID0UT(60) , NCON( 200) , NTERM ( 69 ) , 

X NTRAN(IOO) »NXCOD( 100 ) , POOLED ( 9 ) , REPVAR ( 9 ) , RE SMS ( 9 ) , X ( 99 ) 

DOUBLE PRECISION BO,SUMX,SU MX2 , SUM Y2 , ZEAN 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 


COMMON/SMALL/ 

BYPASS, BZERO 

I, DELETE, FIRST, IFCHI, 

IFSSR, 

12 

X IFTT , 

IFWT, 

INPUT, 

INPUT5, 

INTER, 

13 

X ISTRAT, 

jCOL , 

KONNO, 

LENGTH, 

LIST, 

14 

X NERROR, 

NODEP, 

NOOB, 


NOTERM, 

15 

X NOVAR, 

NPDEG, 

NRES , 

NTRANS , 

NWHERF , 

16 

X P, 

PREDCT , 

REPS, 

RWT, 


17 

X STORY I , 

STORYC, 

STORYX, 

TOTWT , 

WEIGHT, 

18 

X ERRFXD, ECONMY, 

I OUT , ICOL 




19 

LOGICAL ECONMY 





20 

DOUBLE PRECISION RWT , TOT WT , W E I GHT 



21 

LOGICAL BYPASS 

, BZERO, 

DEL t TE 

, IFCHI, 


22 

XIFSSR, IFTT, 

IFWT, 

REPS, 

PREDCT 

, 

23 

xstory:, storyx 

, STORY I 

, FIRST 

, ERRFXD 


24 


COM MON/MAX/M AX PLT 


25 

26 
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c****»»*************************»***»*****************************»»*»»* 


c 


INTEGER CELLS, PLUS 1 

DIMENSION BOUND (451, C ELLBD ( 21 ) , OB S ( 20, 9 ) , RCTI212), 


X STDERR ( 9 ) , VAR(9) , 

X Z(9) 

EQUIVALENCE IVAR, RESMS) 

ODATA BOUND/. 67448907 

1.. 31853932 ,.67448907 

2. . 52439979 ,.84162001 

3. . 43072721 ,.67448922 

4. . 18001235 ,..36610621 
5,1.06756653,1.46521688 

6. . 48877614 ,.67448930 

7. 1.5341 1831.. 1397102 8 

8. . 58945544 ,.76470731 

9. 1.593 23335.. 12566134 

A, .52440000 ,.67448801 
8,1.28154233,1.64490172 / 

LOGICAL SAVRES 
DIMENSION RESPLT ( 1 ) 

DIMENSION SKEW ( 9 ) ,SKUR ( 9 ) 
EQUIVALENCE ( RESPLT, SUMXY) 
DIMENSION I CHAR ( 9 ) , MCHAR (9 ) 
DAT A ( I CHAR ( I ) , I = 1 , 9 ) /6HR ES IDU 
X ,6H VS I , 6HNDEP V,6HAR NO 
DATA ( MCHAR ( I ) , I = 1 , 9 ) / 6HR FS I DU 
X ,6H VS P.6HREDICT ,6HED VAL 
LOGICAL PNCH 


YCALC ( 9 ) , 

YDIFR19) , 

,.43072720 

, .96741604 

,1.15033859 

, .25334708 

,1.28153777 

,.21042836 

, .96741746 

,1.38298075 

, .56594856 

,.79163735 

,.15731067 

,.31863932 

, .88714436 

,1.15034184 

,.28221612 

,.43072722 

, .96741836 

,1.22062834 

,.25334701 

,.38532026 

, .84161868 

,1.03642921 


, 6HALS F0.6HR 
, 1H / 

DEP 

, 6H VAR 

, 1 H 

, 6HAL S FO, 6HR 
, 2HUE / 

DEP 

, 6H VAR 

,1H 


******************************************************** ********** 


JCOL=NOTERM+ NODEP 
NUVAR=N0TERM+1 
BYP ASS = .FALSE. 

KOUNT = 0 
SAVRES=. FALSE. 

I TPLT = NOOB* (NWHERE+2*N0DEP ) 

I F ( I TPLT. LE.MAXPLT ) S AVR ES= . TRUE . 

CALL LRLEGNt I DENT , 54, 0 , . 1, 5 . 0, 0. 0 ) 

CALL LRLEGNI I DENT (10 ) , 24 , 0 , . 1 , 4. 5, 1 . 0 > 


•»*•**••»•»•*•*•*•*••••••••***•*•*•»*••••**•••*•*••*•*••••*•••*»*•**»•« 

I F ( N03 B-2 0 ) 110,120, 120 
110 BYP ASS = . TRUE . 

GO TO 125 
120 CELLS= NOOB/5 

CELLS=MIN0(CELLS,20) 

1= MOO ( CE LL S , 2 ) 

IF(I.NE.O) CELLS=CELLS + 1 
FCELLS= FLOAT ( CELLS ) 

PLUS 1= CELLS + 1 
MINUS1 = CELLS -1 
NDEGCH = CELLS-3 
I R= CELLS/2-1 
IC=IR*( IR-D/2 
IS= IR4-2 
DO 122 J= 1 , I R 
IC=IC*1 
I BC= I S - J 
IRC=IS+J 

CELLBDI IBC)=-BOUND( IC ) 

CELLBO (IRC)* BOUND! I C ) 


27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 
87 
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122 CONTINUE 88 

CELLB3I1) =-l. 0E+37 89 

CELLBO { PLUSl ) =1.0E37 90 

CELLBD ( I S )=0.0 91 

DO 12V K= 1 » NODEP 92 

00 124 1=1, CELLS 93 

OBSII, 0 = 0.0 94 

124 CONTINUE 95 

C 96 

C*»«**««* »*»**■**■»*******«•* 97 

125 DO 133 K=l, NODEP 98 

SKEW (0=0.0 99 

SKUR ( K ) =0 . 0 100 

STOERR I K ) = SORT { ERRMS IK)) 101 

130 CONTI NUE 102 

WRITE(LIST,135) 103 

C 104 

C»***»»***»»*»*»****««*** 105 
DO 433 J= 1 » NOOB 106 

READIINTER) (XI I », 1 = 1,69 ), WEIGHT 107 

IF ( ♦ NOT. S AVR ES ) GO TO 141 108 

INOPLT=NWHERE 109 

DO 143 1 = 1,1 NOPLT HO 

K=( 1-1 )*NOOB+J HI 

140 RESPLT l K) =X( I ) 112 

141 CONTINUE H3 

DO 142 1=1, NOTERM 114 

K= IDDUT ( I ) H5 

X ( I )= XIK) H6 

142 CONTINUE H7 

KBAR=NWHERE H8 

DO 143 1 = 1, NODEP H9 

IC= NOTERM+ I 120 

KBAR=<BAR+1 121 

X ( I C ) = XIKBAR) 122 

143 CONTINUE 123 

C 124 

C****************»***** *************** ********** 125 

DO 163 K= 1, NODEP 126 

YCALCI K ) = BO { K ) 127 

I F ( . NOT. BZERO ) YCALC{K)= 0.0 128 

KBAR= K+NOTERM 129 

DO 153 1=1, NOTERM 130 

YCALCI K J = YCALCIK) + B(I,K)*XII) 131 

150 CONTINUE 132 

ACT0EV= XIKBAR)- YCALCIK) 133 

YDI FRI K ) = ACTDEV 134 

ZIK)=ACTDEV/STOERR(K) 135 

A=ACT0EV**3 136 

SKEWl<)=SKEWIK)+A 137 

SKURK)=SK UR (K)+A* ACTDEV 138 

160 CONTINUE 139 

IFI.NOT.SAVRES) GO TO 179 140 

K=INOPLT*NOOB+J 141 

K8AR = <«-N00B*NQDEP 142 

DO 175 1=1, NODEP 143 

ITC= 1 1 -1 ) *NOOB 144 

I SC = K* I TC 145 

I S = KBA R+ 1 TC 146 

RESPLT I ISC)=YCALCI I ) 147 

RESPLT I IS ) =Z 1 1 ) 148 

175 CONTINJE 149 
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179 CONTINUE 150 

WRITE! LIST, 180) ( X ( K ) , K=NUVAR , JCQL » 151 

WRI TE( LIST, 190) ( YCALC (K ) , K=1 , NODE P ) 152 

WRITE! LIST, 200) 1 YDI FR ( K ) , K=l, NODEP ) 153 

IFIPNCH) PUNCH 5250, J, < YD IFR I K ) , YC ALC I K ) , K=1 , NODEP ) 154 

WRITE1LIST.210) I Z I K ) , K= 1, NODEP ) 155 

IFIBYPASS) GO TO 410 156 

157 

158 

DO 253 K=l, NODEP 159 

DO 233 1=1 »PLUS1 160 

IF ! Z ! < ) -CELLBD! I ) ) 220,220,230 161 

220 OBS!I-l,K)=OBSl I— 1,K)+ 1.0 162 

GO TO 250 163 

230 CONTINUE 164 

250 CONTINUE 165 

166 

410 KOUNT = KOUNT +1 167 

IF(KOJNT.LT.IO) GO TO 430 168 

WRITE!LIST,270) IDENT 169 

KOUNT = 0 170 

430 CONTINUE 171 

******************************************************************* 172 

IF! .NOT.SAVRES) GO TO 439 173 

ITC = NfMERE+NODEP 174 

DO 435 IRC=1, NODEP 175 

IC=!ITC+IRC-1)*N00B+1 176 

DO 434 K=1 ,NWHERE 177 

IS= l K- 1 ) *NOOB+ 1 178 

I CH AR ! 5 ) = I RC 179 

ICHAR! 9 ) =K 180 

CALL L RCNVT ! ICHAR ( 5 ) , 1 , I CHAR 1 5 ) , 1 , 6, 0 ) 181 

CALL L RCNVT! ICHAR !9) ,1,1 CHAR ! 9 ) , 1 , 6,0 ) 182 

CALL LRTLEG! ICHAR, 54) 183 

CALL L RPLOT 1RESPLT! IS ),RESPLTl IC) , NOOB ) 184 

434 CONTINUE 185 

IS=! NrfHERE+IRC-1 )*N00B+1 186 

MCHARI 5 ) = I RC 187 

CALL L RCNVT t MCHAR l 5 J , l.MCHAR! 5) , 1,6, 0) 188 

CALL LRTLEG(MCHAR,54) 189 

CALL LRPLOTIRESPLTIIS) ,RESPLT( IC) »NOOB) 190 

435 CONTINUE l g l 

439 CONTINUE l g 2 

193 

t****. ******** *********** ********************************************** 194 

IFIBYPASS) RETURN l g 5 

DO 653 K=l, NODEP 196 

SKEW(<»=SKEW!K)**2A! FLOAT! NOOB ) **2*ERRMS < K ) **3 ) 197 

SKUR!<)=SKUR«K)/! FLOAT l NOOB ) *ERRMS I K ) **2 ) 198 

CHI S0= 0. 0 199 

DO 643 1=1, CELLS 200 

RCT ! I ) =QBS ! I , K ) 201 

CHISQ=CHISQ+RCT!I)**2 202 

640 CONTINUE 203 

CHISQ=FCELLS*CHISQ/FLOAT 1 NOOB ) -FLO AT ! NOOB ) 204 

WRITE! LIST, 280) NDEGCH , CHI SQ , SK EW < K ) , SKUR ( K ) 205 

CALL HIST!K, RCT, CELLS) 206 

650 CONTINUE 207 

RETURN 208 

135 FORMAT! 51H FOR EACH DEPENDENT TERM AND OBSERVATION IS PRINTED 209 

X 7 31H OBSERVED RESPONSE !Y OBSERVED) 210 

X 7 29H CALCULATED RESPONSE IY CALC) 211 
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X / 28H RESIDUAL IY08S- YCALC=YDIF ) 212 

X / 28H STUDENT I ZED RESIDUAL ( Z ) ) 213 

180 FORMAT ( 12HKY OBSERVED .9G13.4) 214 

190 FORMAT ( 12H Y CALC ,9G13.4) 215 

200 FORMAT ( 12H Y DIF f 9G13.4) 216 

5250 FORMAT (I6»4E16.8/(6X»4E16.8) ) 217 

210 FORMAT ( 12H STUDENTIZED .9G13.4) 218 

270 FORMAT (1H113A6,A2) 219 

280 FORMAT { 27H1CH I -SQUAR E STATISTIC WITH I6.22H DEGREES OF FREEDOM = 220 

X GL4.6/12H SKEWNESS = G14.6/12H KURTOSIS = G14.6I 221 

END 222 


CRSPLT PROGRAM 


CRSPLT accepts a subset of the data used for a NEWRAP problem. It can be used 
as a preregression analysis program to help formulate model equations to be analyzed 
with NEWRAP or it may be used as a postregression program by using punched output 
from NEWRAP to obtain more complex residual plots than direct use of NEWRAP allows. 

When used as a preregression program, it will compute an X'X and C matrix in- 
cluding all the terms (independent and dependent) if requested. It can also compute ei- 
genvalues and eigenvectors of the submatrix of X'X corresponding to the independent 
variables. 

When used as a postregression program, the punched output of residuals and pre- 
dicted values from NEWRAP can be plotted against new functions of the independent 
variables. 

The input is much the same as for NEWRAP. The seven sets of input are as fol- 
lows: 

(1) IDENTIFICATION (I, IDENT)(I2, 13A6): EDENT is Hollerith data used to iden- 

tify the problem. I indicates the number of additional cards to be read for iden- 
tification (columns 1 to 78). 

(2) PROBLEM SIZE (NOVAR, NODEP, NOTERM, NOOB)(3I4, 15) 


NOVAR 

NODEP 

NOTERM 

NOOB 


Number of input independent variables 
Number of input dependent variables 
Number of terms in model equation 
Number of observations 


(3) TRANSFORMATION S: This input is the same as the transformations of NEW- 

RAP except for the upper limit of 150 transformations. 

(4) FORMAT (INPUT, FMT)(I2, 13A6): INPUT specifies the unit number the input 

data is stored on and FMT indicates the format for reading it. 
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(5) PLOTTING REQUESTS 

NOPLTS (14) 

(IXPLT, IYPLT) 

One card supplies NOPLTS, the number of plots desired. The following 
cards supply pairs of integers indicating which pairs of terms to plot. The for- 
mat is 4012 (i. e. , 20 plots per card). The first number of the pair (IXPLT) 
specifies the sequence number of the term to be used as the abscissa. The 
second number (IYPLT) specifies the sequence number of the term to be used 

as the ordinate. As an example, the following sequence of transformations and 

2 

subsequent plotting requests would cause X-, to be plotted against x., as well 

3 *■ 

as against x^: 

MODEL SIZE 0003 

TERMS 616263 

TRANSFORMATIONS 01000061 01026162 01026263 

CONSTANTS blank card 

NOPLTS 0002 

(IXPLT, IYPLT) 0102 0302 

(6) M ATRIX REQUESTS (XTXC, EIGENC)(2L1): If XTXC is F, no matrix calcula- 

tions are executed. If it is T, than an X’X, X'X deviation and a correlation 
matrix of all the NOTERM + NODEP terms appearing on the TERMS card are 
computed. If EIGENC is T, then the eigenvalues and eigenvectors of the sub- 
matrix corresponding to the independent terms (the first NOTERM terms) are 
calculated. 

(7) DATA: Same usage as in NEWRAP. 

An illustrative set of input is given, followed by the corresponding sample of output 
and a main program listing. The subprograms TRIANG, RECT, and EIGEN are re- 
quired and are the same as in NEWRAP. 

1 *> SAMPLE CRSFLT PROBLEM 

DATA IS FROM DRAPER AND SMITH 

APPLIED REGRESSION ANALYSIS (REFERENCE 4 OF NCHRAP REPUPTi 
CHAPTER 7 

INITIAL MODFL Eu (sAS 

Y = CHAMBER PRESSURE 

X 1 = TEMPERATURE OF CYCLE 
X2 = VIBRATIUN LEVEL 
X3 = DR0P< SHOCK) 

X4= STATIC FIRE 

Y = BC <■ B 1 X 1 <- B 2X2 + B3X3 + D4X4 «■ ERR 
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THE FOLLOWING TERMS ARE BEING CREATED FOR RESIDUAL PLOTS 

1 XI 2 X2 3 X3 4 X4 

5 X1*X2 6 X1*X3 7 X1*X4 8 X2*X3 

4 X2*X4 1(5 X3*X4 11 YIPRED) 12 UNIT NO. 

13 RESIDUALS 

6 1 12 24 

1 3 

6 162 £3 646 56 6 6 7686970 717273 

'''imnf!72T’20r0f!61030r0C62f 40C tC6 305QG0O 6406Q0C073' , )7QC ! rP7 1010 2626501? 26366? 1? 2646 7 

>' 2C 2636 80 20264690 30 26470 
03 (8X.F2.0,2X.4F6.0/6X.2E16.8) 

12 

01 1 3C21 303 1 30413051306130 7 130 81 3091 31C 13 11 131 21 3 
TF 

UNIT NO. I -75 0 0 -65 1.4 

1 -n. 3655C47OE+01 C . 50 5 50470E +C 1 

2 175 0 C 150 26.3 

2 -0. 76331902E+C0 0. 270 63 31 9E +02 

7 ( 0 -65 15C 29.4 

3 O.23366811E+01 f . 27063 31 9E +02 

8 0 C 165 -65 9.7 

4 0.464 49 530E+C1 0 . 505 50470E +01 

9 000 1 50 32.9 

5 0. 58366811E+C1 f . 2 70 6331 9E +02 

10 -75 -75 0 150 26.4 

6 O. 50664449E+CC C . 2 5893355E +C2 

11 175 175 C -65 8.4 

7 0.61503899F+00 0 . 7 7849 61QE +C 1 

14 -75 -75 -65 15C 28.4 

8 0.25066445E+C1 C. 25 893355E +02 

15 175 175 165 -65 11.5 

9 0.3715039CE+C 1 O . 7 7849 610E +C1 

18 o 0 -65 -65 1.3 

10 -C .37550470E+O1 f . 50 55G 47CE +C 1 

19 O C 165 150 21.4 

11 -0. 56.6331 89E + 01 0 . 2 706 331 9E +0 2 

20 0! -75 -65 -65 0.4 

12 -0.348 50 838F+C1 0 . 3 88 50 83 8E +0 1 

21 0 175 165 15C 22.9 

13 -O.68932328E+01 (. . 29 793233 E+C 2 

24 non -65 3.7 

14 1355O470E + O1 C . 50 5 50 4 70E +0 1 

3 0 -75 C 150 26.5 

15 L . 6t~ 6 644 39E +0C C . 25 893 355 E +0 2 

5 O -75 O 150 23.4 

16 -0.24933555E+C1 0 . 25893355E +02 

16 o -75 0 15C 26.5 

17 0.606 644 39E +00 f . 2 5 893 355 E +C2 

4 Oi 175 O -fc5 5.8 

18 -f. 1984961OE+0 1 C .77849610E+C1 

6 Oi 175 r -65 7.4 

19 -0. 384961C1E+0C 0 . 7 7849 610 E +C' 1 

17 0 175 O -65 5.8 

20 -O. 19849610E + 01 C . 7 784961'~E +C 1 

12 O -75 -65 - 1 5C 28.8 

21 0.29C66443E+01 C . 2 5 893 355 E +C2 

22 o> -75 -65 -150 26.4 

22 Q • 50 6 6444 9F + t C 0 . 2 5 893355 E +02 

13 ^ 175 165 -65 11.8 

23 0.401 50389E+01 0 . 7 784961C E +0 1 

23 o 175 165 -65 11.4 
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1 1 R F T C CRSPL X 


CGMM0N/B1/ X (99 I ,CUN( 991 , SUMX ( 70) , SUMXXI2485 ) . AI70 ,70) 1 

X. XI) AT A ( 1 20 rtf' ) 2 

COMMON /B2/XME AN (701 ,XSTD( 7(. » . SUMX 2 (70, NT RAN S , NCON ( 3^0 ) , 3 

X NTERm(70J , NTRAN ( 1 5< ) , NXCUO ( 1 50 I 4 

0 1 MENS ION IDENTI 13) ,FMT( 13 l.FMTTR I ( 141 ,CORR( 1) .FMTSGH3) 5 

QATAtFMTTR I < I) . 1 = 1, 3) /6HI 1H I 6 , 6H , 8G1 5 . , 2 H6 I / 6 

DATAt F mTSGL ( I) . 1=1, 31 /6H< 1H I 6 , 6H , 8G1 5 • , 2H6 ) / 7 

LOGICAL XTXC.EIGENC 8 

EQUIVALENCE (CORK. A) 9 

DATA (CHAR ( 2 I/6H VS / 13 

COMMON IXPLT (4001 , IYPLT(4CC ) 11 

DIMENSION I CHAR ( 3 ) 12 

C ****** ****** ******************** ************ *************************** 13 

lO READ (5.110) I, IDENT 14 

wKITE(fc.lU) I OENT 15 

DO ICO J=l, 19863 16 

X ( J 1=0.0 17 

1 Of* CONTINUE 18 

113 l F( I) 120. 120,115 19 

115 K FAO( 5.300 ) EMT 2C 

WPITE(fe.JOi) fmT 21 

1 = 1-1 22 

GO TO 113 23 

120 RFAD(5,112) NO V AR . NUDE P . NO TER M , N JOB 24 

*R I TF< 6. 30 5 ) NO VAR .NODE P. NUT ERM.NUOB 25 

MVTERM=NOTERM+NODEP 26 

L=MvTEKM*N()OB 27 

1 Ft L.GT. 12000 GO TO 10C0 28 

RFAU (5.282) NT RANS , KCNNU 29 

I F ( NTR ANS ) 2 5 5 . 2 55 . 2 2C 33 

220 REAG( 5,230) (NTERM( K ) ,K=1 .MVTERM) 31 

w R ITE(6,235) (NTERM( N ) ,K = 1 , MVTERM) 32 

REAO( 5,2 30) ( NX COD ( I) . NTRAN ( I ) , NCOM 2*1-1) , NCON (2* I I , 1 = 1 , NT RANS ) 33 

Kk (TE( t>. 240 ) (NXCUO ( I ) , NTRAN ( I ) ,NC0N(2*I-1) ,NCON< 2* I ) , I =1 , NT RANS ) 34 

IF(KCNNO) 2 55.255.25C 35 

2 50 R t AD ( 5,2 60 ) { CON (It, 1 = 1, KUNNO ) 36 

WR I TF< 6 , 2 62 I (CCN ( I I , I=1,K0NN0) 37 

2 55 k EA D ( 5 . 3000 I (UIN.FMT 38 

w R I T E ( 6 . 30 0 1 ) IUIN.FMT 39 

RFAD( 5,282 ) NOPLTS 43 

I F ( NOP LTS.LE.O) GO TO 320 41 

IF ( NOPLTS. GT.3C0) GO TO 2000 42 

K F A D< 5 , 2 30 ) ( IXPLTI I ) , I YPLT ( I), 1=1, NOPLTS) 43 

wR ITF( 6. 5000 ) ( IXPLT(I),IYFLT( l), 1=1, NOPLTS) 44 

320 CONTINUE 45 

RFAD( 5,4000 ) XTXC.EIGENC 46 

NDRD=NOVAR*NCDEP 47 

C 48 

C ***************** **************************************************** 49 

RNOOB = l.n /FLOAT ( NOOB ) 50 

00 690 J = 1 , NOOB 51 

READ( IUIN.FMT) ( X ( l ) . I =1 . NURD ) 52 

I F( NTRANS) 450,450. 34C 53 

340 CALL TRANS 54 

00 430 K=l, MVTERM 55 

1= NTFRM(K) 56 

X ( K ) = C 0 N ( I ) 57 

430 CONTINUE 58 
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4 5 r CONTINUE 59 

IF! .NUT.XTXC) GO TO 620 6P 

l R=0 61 

DO 610 1=1, MVTERM 62 

SUMXIII =SUMX( I I + X(I) 63 

DO 600 11=1,1 64 

I R = IR + 1 65 

SUMXX 1 IR I = SUMXXIIR) * X(I)*X!II) 66 

600 CONTINUE 67 

610 CONTINUE 68 

620 CONTINUE 69 

I A = J— NOOB 7C 

DO 650 I = 1 , M VTE RM 71 

IX = IA ♦ l*N003 72 

X DATA I I X ) = X ( II 73 

650 CONTINUE 74 

690 CONTINUE 75 

I FI .NUT.XTXC ) GO TO 72C 76 

WRITFI6. 101G ) 77 

CALL RECTINOOB. MVTERM. NOOB, MVTERM, XDATA.FMTSGL) 78 

WRITE! 6. 1020 ) 79 

CALL TR IANGISUMXX, MVTERM. 8, FMTTRI ) 80 

I R= 0 81 

00 700 1=1, MVTERM 82 

IR = IR + I 83 

SUMX2! I»= SUMXX ( I R I -SUMXI I )**2*RN00B 84 

700 XME AN 111= SUMXI I )*RNOOB 85 

I R = 1 86 

DO 710 J=l. MVTERM 87 

DO 710 K=1»J 88 

SUMXX1 IR I = SUMXXIIR) - SUMX I K ) *SUMX I J ) *RNOOB 89 

CORRIIR) = SUMXXIIR)/ SQR T 1 SUMX2 I J ) * SUMX2 I K ) > 90 

I R= IR + 1 91 

710 CONTINUE 92 

WRITE16, 1025 ) 93 

CALL TR I ANGISUMXX, MVTERM, 8, FMTTRI ) 94 

WRITEI6.1030 ) 95 

CALL TRIANGI CORR, MVTERM, 8, FMTTRI ) 96 

IFI .NOT.EI GENC ) GO TO 720 97 

CALL EIGEN! SUMXX, A. NOTERM, C ) 98 

WR ITEI 6, 1040 ) 99 

J = 0 ICO 

DO 718 1=1. NOTERM 191 

J= J* I 192 

7 1 R SUMXX! I ) =SUMXX I J) 1C 3 

WR ITEI 6, 1050 M SUMXX! I ), I =1 .NOTERM) 1C4 

WRITE16.1060 ) 195 

CALL RECTI NOTERM, NOTERM, NOTERM, NOTERM, A, FMTTRI ) 1C6 

720 CONTINUE 197 

CALL LRLEGNI IDENT, 54,0,. 1 ,5.0, C.O) 1C 8 

CALL LRLEGNI IDENTI 1C), 24,0, .1,4. 5, 1.0 199 

DO 900 K= 1 , NOPL T S 119 

I CHAR I 1)=IXPLT I K) 111 

ICHAR(3)=IYPLTIK) 112 

IS1=( IXPLTIK )-l )*N00B + 1 113 

IS2 = < IYPLTI K )-l )*N0QB + 1 114 

CALL LRCNVT I ICH AR II), l.ICHARIl), 1,6,0) 115 

CALL LRCNVT ( ICH AR (3) ,1, ICHARI3 ), 1,6,0) 116 

CALL LRTLEG! ICHAR.18) 117 

CALL LRPLOT I XDATAI IS1I.XDATAIIS2), NOOB ) 118 

900 CONTINUE 119 

GO TO 10 12 r 
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If ■)(' wk [ T F 1 B • 95 ) L 121 

STOP 122 

9 5 FOR MAT 1 62H T Ft REQUIRED NUMBER OF LOCATIUNS EXCEEDS THE 1200" AVAI 123 

XLABLE 18) 124 

?r-rr wr j TE< 6. 2fr 5 > NOPLTS 125 

STOP 126 

?! t 5 FORMAT ( 2 5H MAX NO. OF PLOTS IS 3C0 18) 127 

IIP FORMAT ( ( 2. 13A6 I 128 

111 FORMAT ( 1H1 . 13A6 ) 129 

1 1 2 FORMAT ( 3 14. IE) 13C 

4"uP FORMATI2L1) 131 

5GCC FORMAT! 1HK/ I 151 IX. 12 . IX. I 2 ,2X) ) ) 132 

3UO FORMAT ( 1 3A6, A2 ) 123 

301 FORMAT (1H 13A6.A2) 134 

305 FOKMATIlH 3161 135 

282 FORMAT! 214) 136 

2 30 FORMAT! 40 12 ) 137 

236 FORMAT < 1 1H TFRMS ARE / !1H 30141) 138 

241 FORMAT < 25H TFE TRANSFORMATIONS ARE /!1H 51414, 5X)») 139 

26C FORMAT! 5E15.7) 14" 

2b? FORMAT I 19H THE CONSTANTS ARE / t ! 1H 8G15.7))) 141 

1C 10 FORMAT! 16H1 TFE DATA MATRIX ) 142 

If 20 FORMAT! 26H2TFE X TRANSPOSE X MATRIX ) 143 

1325 FORMAT! 32H X TRANSPOSE X DEVIATIONS MATRIX ) 144 

1030 FOkMAT! 23H2TFE CORRELATION MATRIX ) 145 

If, 40 FORMAT ( 46H2 F CL L Ort I NG ARE EIGENVALUES OF X TRANS X MATRIX) 146 

1050 FORMAT < 1 H 8G16.7) , 147 

106r FORMAT! 53HKF IGENVECTORS BY COLUMNS IN SAME ORDER AS EIGENVALUES) 148 

30L.fi FORMAT! 12. 13A6) 149 

3001. FORMAT 1 1 H I6.2X.13A6) 15" 

F NO 151 


11RFTC TRANSX 

SUBROUTINE TRANS 1 

O********^******** *$>!=** ****’(‘*******4************** ******* ******* ******** 2 

C 3 

COMMON /H 1 / X 199 ) .CON! 99) , SUMX! 70) , SUMXXI2485 ) , A< 70,70) 4 

X. XDAT A ( 12000 ) 5 

C0MMGN/B2 /XMEAN170) , X S TD ( 7C ) , S UMX2 1 70 ) , NTRANS, NCON I 3"0 ) , 6 

X NTFRM170 I .NTRAN ( 1 5C ) , NXCOD! 1 50 ) 7 

C 8 

f. **** **** * * *** **************** *** ** * ** *** **** ** ** * ** ************ ***** *** g 

C THIS SUBROUTINE PERFORMS TRANSFORMATIONS IF THIS OPTION IS 1" 

C REUU6STEC. 11 

C. 12 

C 13 

f. K TRANSFORMATION SET NUMBER. 14 

C NCON ( 2 *K— 1 ) CONSTANT NUMBER TU USE. 15 

C N CUN ( 2 *K ) DERIVED CONSTANT. 16 

C. NTRAN! M NUMBER OF TRANSFORMATION REQUESTED. 17 

C NXCOOIK.) VARIABLE NUMBER 18 

r. 19 

8C DO 5"0 K = 1 , N TRANS 2" 

I =N CON ( 2 *K- 1 ) 21 

l FI I I IPO . IOC .IIP 22 
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l r ‘ con 5 = ‘.i ?3 

Go TO 120 24 

11 ” cow s=con < I ) 25 

l?” l = i\iXCOI)< K ) 26 

Y = X ( I ) 27 

MT8AN = NTRAMK) 28 

I F( MTRAN.Lt.” I MTKAN=32 29 

14( GO TO< 150,160,1 7P,16C,19r,2C” ,21” ,22(.',23”,24P,25”,26''',27: ,28C,29”, 3” 

X3 r ” ,310, 32” , 33” ,34” ,3 50 ,36” ,37C ,3 8<' ,3 9” ,440,41” ,42”, 43” ,44” , 31 

X 442.45' ) , M TRAM 32 

1 5f CUN S= Y +CUNS 33 

GO TO 46 ,f 34 

1 6< C GNS = Y 4C0N S 35 

GO TlJ 460 36 

17. CUNS=CONS/Y 37 

GO TO 460 38 

18” CON S= FXP ( Y ) 39 

GO TO 460 43 

190 CONS=Y**CONS 41 

GO TO 460 42 

2 0t. CGNS= AlflGl Y ) 43 

GO TO 460 44 

21t C0NS = ALCJGlO( Y) 45 

GO TO 4o0 46 

??(■ C 0 N S = S I N ( Y I 47 

GO TO 460 48 

2 3” crjNS = cns<Yi 49 

GO TO 460 5” 

24” C(JNS = SIN< 3. 141 55265*<CCNS*Y) ) 51 

GO TO 460 52 

2 50 r.ONS=r.US< 3.14159265*(CGNS*Y) ) 53 

GO TO 460 54 

26 r 0005=1.0/7 55 

GO TO 460 56 

27 r CUNS=EXP(CGNS/Y > 57 

GO TO 460 58 

28f CONS=FXP(CCNS/(Y*Y) I 59 

GO TO 460 62 

29” CuNS=SURT< Y) 61 

GO TO 460 62 

3 OF C(/NS=1.T / S (J R T ( Y I 63 

GO TO 460 64 

3 If CnNS=OONS**Y t 5 

GO TO 460 66 

320 C0NS = 1”.”<**Y 67 

GO TO 460 68 

3 3' C 0 N S = 5 [ N H ( Y I 69 

GO TO 46 0 7” 

3 40 CONS=COSH(Y) 71 

GO TO 4fc>n 72 

350 CONS = < l.o-CD5( Y) )/2.0 73 

GO TO 460 74 

3 60 CONS=ATAN< Y ) 75 

GO TO 4t>0 76 

370 CUN S = A T AN 2 ( Y /CONS ) 77 

GO TO 460 T 8 

38” CONS=Y*Y 79 

GO TO 46” 8” 

391 CUNS=Y*Y*Y 81 

GO TO 46” 1 82 

4 CUN S = AkS IN < SURT ( Y I 1 83 

GO TO 4fcO 64 
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I 


t r. 

CONS=2.0*3. I *159265*Y 

85 


GO TO 460 

86 

A 20 

r.GNS = L.O/( ?.C*3.1A15S265*Y) 

87 


GO TO A60 

88 

4 3C 

CON S = ERF ( Y) 

89 


GO TO 460 

9r 

440 

CONS- GAMMA ( Y I 

C 1_ 


GO TO 460 

92 

AA? 

C£«NS = Y/CONS 

9? 


GO TO 460 

94 

A 50 

CON S*Y 

9? 

A 60 

I *NC0N(2*K ) 

°6 


I F f I )47Q , 470. 4B0 

97 

A7F 

CONI K ) =CON S 

°8 


GO TO 500 

99 

ABC 

CON ( I ) =CQNS 

KC 


I F ( 1-60) 500.500,490 

ir 

49C 

XII) ■CONS 

I r 2 

500 

CONT INUE 

103 


RETURN 

l r A 


END 

1:5 


SAMPLE CftSPLT PROBLEM 

DATA IS PROM DRAPER AND SMITH 

A PPL I EO REGRESSION ANALYSIS (REFERENCE A OF NEWRAP REPORT) 
CHAPTER 7 

INITIAL MODEL EU wAS 

V= CHAMBER PRESSURE 
X l a TEMPERATURE OF CYCLE 
X 2= VtHRAT ION LEVEL 
X3 = CROP ( SHOCK) 

X4 = STATIC FIRE 

Y = HO + BlXl ♦ tt2X2 B3X3 + B4X4 ♦ FRR 
THE FOLLOWING TERMS ARE BEING CREATED FUR RESIDUAL PLOTS 


1 XI 

2 

X2 


3 

X3 



A 

XA 




5 X 1*X2 

6 

X1*X3 


7 

X 1 *X A 



8 

X2*X3 




9 X2*XA 

1C 

X3*XA 


11 

Y ( PKED ) 


12 

UNIT NO. 




13 RES IOUAI S 












6 

l 12 












? A 













TERMS ARF 













61 62 

63 6 A 65 

66 67 

68 69 

7C 

71 72 

73 







THE TRANSFORMATIONS ARE 











1 C 

0 72 

2 C 

C 61 


3 

o r 

62 


A '' 

r 63 

5 r r 

6A 

6 r 

O 73 

7 0 

C 71 


1 

2 62 

65 


1 2 

63 66 

1 2 6A 

67 

2 2 

63 68 

2 2 

6A 6 9 


3 

2 6A 

70 






5 

I QX ,F2.0. 2X . 

4F6.0/6X, 

.2E16.8) 










1 13 

2 13 3 13 

A 13 

5 13 


6 13 

7 13 

8 

13 

9 13 

l r 13 11 13 

12 13 



THF DATA 
1 

MATRIX 
-75. Of 000 

C 

0 

-65. 0° 0C? 

c 

n 



2 

175 .0-000 

r 

c 

15 r .CC? r> 

c 

r> 

-a y n r- r r 


3 

0 

n 

-65. 

1 5 r . n pr c 

r 

-A55.oom 

I “5 "."or 


A 

c 

r 

165. •'OCC 

-65.CCPCC 

r 

13?o.rrr 

_* ? r # rrrr 

" 

5 

0 

r 

0 

150.C07P 

C 

r* 

1 ■» 6". 

‘ 

6 

-75.00000 

-75. room 

0 

iso. 0 °-:? 

-75?. or : ~ 

r 

! 50 *- . - 0 - 

- r 

7 

175.0000 

175.0000 

n 

-65. P n O r ' n 

192 5. r P‘? 

r» 

-715. 

~ 

R 

-75.orooo 

-75. ooroo 

-65.0COOO 

15C. p C20 

- 1050 . rc^ 

-91-.-O0- 

21-''', o'-- 

65 7 r , - - - 

9 

175.0000 

175.C0OC 

l65.0 n 0r 

-65.P n P3? 

2625. TO" 

2A75.T7- 

-975. " r ’* 

, ’88" r ' ; . 

10 

o 

C 

-65.0PO0C 

-65.orn'0 

r 

- 1 1 7 r • 0-0 

- 1 1 7°. ''nr 


11 

0 

0 

16 5. OP 00 

15P. r C r ;0 

r 

31 3 5. 0 0 r 

2 P 5" .-o'* 

■* 

12 

0 

-75.ccrce 

-65.PCP3C 

-65. orn-* 

-15"". r^.' 

- 13 - 0 . 0 -'' 

-l 1 "*."'' 


13 

r 

175.0000 

165.30OC 

lSC.^rr 

3675.0 C" 

346*. ro-> 

315". o- 

" 

1A 

O 

0 

o 

-65 .PoC?" 

n 

r 

-156".-"" 

" 

15 

0 

-75.00000 

0 

i5 rt . n c:p 

-225. COP" 

0 

A 5-. " -o “ 

" 

16 

0 

-75.0LCGC 

p 

15P. r P9r 

-375. Cf r r 

r 

7 C -. --or 


17 

n 

-75.corcc 

r 

1 5P. r DPC 

-120r.( nr 

r 

26''-. or - 

“ 

18 

0 

1 75.CCOC 

n 

-65.o n o:c 

7C9. rrrr 

*" 

- 28 "'.o**o.* 

' 

19 

0 

175.CO0C 

0 

-65 D M ? 

i 05 r.ro- 

r 


- 

20 

o 

175.OCO0 

0 

-65 .r^c^c 

2975. Co? 

0 

-! nc . """ 

* 

21 

o 

-75. cor 00 

-65.o^eoc 

-15 O.rr^n 

- 90 ?, 

-79 n . ooo r 

-18"'".""" 


22 

0 

-75.OOC00 

-65.0Pnor 

-i5 /, .pc , rc 

-1650. rr 

- 14 30 . rrr 

- 3 ? o'* • - o- 

- ~ 

23 

0 

175. omr 

165.CP0C 

-65.r p ??? 

2275. oor 

21A5. rr o 

- 8 A 5 • - o " 


2A 

r 

175.000C 

165.0CO'' 

-65.oc n :c 

4?25.C r 

37 7 5 , 0 “ 0 

_i 495 . 00 ' 


1 

2 

3 

A 

5 

6 
7 
6 

A875.T OO 
26250.00 
0 

-0 

0 

-11250.00 
-11 375.00 
-11250. CO 

-0 

c 

r 

-0 

r 

-11250. Cn 
- 1 1 3 75. n 0 
-1125C.OC 

5.C550A7 
27.06332 
27.06332 
5.C55C A7 
27.06332 
25. 89335 
7. 78A961 
25.89335 

2m r 'f'C0 r C 

7. error? 

8. r>rr^f >7 

9. poc ?r c 

10 . rpcPC 

11. c''’'': 
14 .rr^rr 

-3.655047 
-C.763319 
2.336681 
4.64495? 
5. 836681 
C .576644 
r.615039 
2.5 "66^4 
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9 

-11375.00 

-11375.00 

7.784961 

i5.o?^:o 

3.715 r 39 

lo 

-r 

-r 

5.C55°47 

18. po??? 

-3.755 r 47 

11 

r 

c 

27.06332 

19.C n 3?0 

-5.663319 

12 

-0 

4875. PCC 

3.885C 84 

2 C. 0 ??v? 

—3. 48 50 8 ^ 

13 

o 

26250.00 

29.79323 

2i. COPCO 

-6.893233 

14 

-( 

-C 

5.055047 

24.P0C0C 

-1.355347 

1 5 

r 

-11250.00 

25.89335 

3.CC000O 

C. 676644 

16 

n 

-11250.00 

25.89335 

5. OC^OOC 

-2.493356 

17 

n 

-11250.00 

25.89335 

16.0030" 

0.606644 

18 

-0 

-11375.00 

7.784961 

4.orroco 

-1.984961 

19 

-0 

-11375.00 

7.784961 

6.007330 

-C.3B4961 

20 

-r 

-11375.00 

7.784961 

17.0003'' 

-1.984961 

21 

-n 

1125C.r.O 

25.89335 

12 . 0000 ? 

2.906644 

22 

-o 

1 12 5C. or 

25.89335 

22 .COOCC 

C. 576644 

23 

-o 

-11375. ''C 

7.784961 

i3.oro:o 

4*01 5 A 39 

24 

-0 

-11375.00 

7.784961 

23.CO0CC 

3. 61 5039 


X 

TRANSPOSE X MATRIX 







1 

108750.0 








2 

72500.00 

290 OCC. 0 







3 

33750.00 

13 5000. C 

1887CO.C 






4 

-14125.00 

-82250.00 

15050. OC 

320700.0 





5 

931250.0 

C.39425CE+C7 

0.2410 5CE <-0 7 

-521125. C 

0.6 50C CO E+? 0 




6 

501375.0 

C.241C5PE+07 

C.30882CE+07 

644525.0 

0 • 4 5 352 5E+C 8 

A . 554528E+"8 



7 

-50H375.0 

-521 12 5. C 

644525. f 

0.3826CCE+C7 

- 0 . 175562E+C 7 

0 . 16 r 57 1 E +" 8 

r . 5Q9Q37F+ r 8 


8 

0. 468750E *07 

0.46875CE+07 

0.44475CE+C7 

-C.114562E+07 

C.7C6781 5 + 0 8 

0 . 67 r 294E +" 8 

- r . 179156E+ A P 

a575’lF + ’ , o 

9 

0. 1 93437E +07 

-0.229375E+07 

-C • 114562 E+0 7 

0. 172438E+P7 

-C. 3 15 r 62E + T 0 

- A . 1 79 15 6 E + r 0 

- r . 137! PlE+^fl 

- r . 3 P3297F + ~q 

10 

-0.2293 75E+07 

-C.717500E+C7 

-0.23475CE+07 

-0. 30 1625E*r 7 

-C.761969E+C 8 

-O.257269E+O0 

3 7 7037E+ r 7 

? 8?2°7E+~9 

11 

3197. 685 

865. T004 

6679.541 

27249. 98 

65669.14 

156194. 1 

271582.'' 

?*» 

12 

3025. non 

11600. OC 

10290. OC 

4CC.OOOC 

1 98 5PC • 0 

191380.'* 

-4465 r .?••' 

cm 17 

13 

672.3147 

-C. 1 52 58 8E-C4 

500.4591 

-1023.986 

423.3643 

- 1394. r 72 

-2991''. 86 

1°4'5' .6 


9 

0.122473F+1P 








10 

0.511 9C6E+C9 

0.250450E+10 







1 1 

-246 52.88 

-692766. 1 

9052.137 






12 

-508375.0 

-521125.0 

4778.626 

49?0. OCC 





13 

-121009.6 

-265546.4 

0.6 86646E-04 

-82.22644 

256.3119 




RANSPQSE X DEVIATIONS MATRIX 







1 

105000.0 








2 

62500.00 

263333.3 







3 

26250.00 

115000.0 

173700. ? 






4 

-19375.00 

-96250. OC 

4550.000 

313350.0 





5 

786250.0 

0. 3555 83E + C7 

C. 21205CE+07 

-724125.0 

0.593933E+ A 8 




6 

372750.0 

0. 20675CE+07 

0.2 83C95E+0 7 

46445C.C 

C.4?379 r E+08 

0. 51 r 91 r E +0 8 



7 

-513375.0 

-534458.3 

634525.0 

0. 3819CCE+C 7 

-0.194896E+C 7 

A . 158856E+ r '8 

C . 5°9P71E+ r P 


a 

0.4265626+07 

0.35625OE+C7 

C.360375E+C7 

— C. 173625E+0 7 

C‘.593656E+ r 8 

C. 52 5591E +0 8 

-?. 189781E+'*8 

". Bl‘*' , 7"F + "0 

9 

0.21ir94E*r7 

-C.182292E+C7 

-7925C0.0 

C. 197156E+0 7 

-0. 246792 E +08 

-0. 118595E+0 8 

- r . 1348278+^8 

363439F + ** ° 

10 

-0.1265636+07 

-0.443333E+C7 

-2 91250.0 

- A .157688b+C7 

-C. 3 64927 E + 0 8 

?• 953781E + r 7 

r . 5!5 r 21E+ A 7 

?676”F+"9 

1 1 

-1776.064 

-12398. 33 

-3267. 958 

20286.74 

-126649.2 

-14405.53 

22495'*. 9 

-2"“8 526. r 

12 

-725.0000 

160C.CC0 

279C.C0P 

-485?. OCC 

5350C.P r 

62755. ? A 

-4Q65'* » rr 


13 

672.3146 

-0. 140429E-03 

5CP.4590 

-1023.986 

423.3624 

-1394.073 

-2 Q 91 ” . 56 

110491 # A 


9 

10 
1 1 
12 

P. 121642E+10 
0.463499E *09 
209527. 8 
-331812. 5 

r.222 262E+10 
67C87P.2 
51*7000, P 

2455.289 

-195.1232 

115 A . Of? 


13 

-121009.6 

-265546.3 

C.640825E-C5 

-82.22649 

256.3119 


f.DRRFLAT ION MATRIX 

1 l.ncroco 

? 0.375865 

3 0.194372 

4 -0.106815 

5 0.314845 

6 0. 161014 

7 -0.2C4556 

l.roocm 
C. 537706 
-0.335C 68 
0.899124 
0. 563940 
-P.134472 

l.PCOCQC 
0.19 5 r '28E-Cl 
C. 660191 
C. 950764 
C. 196571 

1 . CONOCO 
-?. 167853 
C. 116135 
0.880858 

n Aprn«“ 
0.733377 
-?.326516E- A 1 

1.0P r P? r 
C • 2 87 r P 8 

1_ , 

■»** “*r <-r 


8 

0.462516 

P.243916 

P. 3 A 2804 

-r. 8977 

0.297853 

r . 258^8 A 


tp q ? ? qc-* - ^ 

1 . 

9 

0.186784 

-f • 1C 1 8 53 

-C. 5452026-01 

r . 160984 

— C .91 8163E-0 1 

- r . A 759 ? 5 F-''? 

_r e 

AO Qi_ 2 1 c - r ’ 

"> A91 1 Q 

10 

-0. 8289716-01 

-0 • 1 83250 

-P. 148229E-01 

-■?. 597516E- r 1 

-?• 1 '*'* 302 

r . 2 p 3176E-''l 

' . 

191 

- *. I ooASf 

11 

-0.110615 

-0.487595 

-C. 158243 

'*.731385 

-C. 331652 

Uf 9^2 QC_r> i 

" . 

58M4** 

-' . 1 47P^O 

12 

-0.6597736-01 

0.9 19429E-P 1 

0.1974*4 

- A . 255492 

0 • 2 r 97? 9 

?. 259F29 

-' . 

1 av** 2 C 

- o P 2 j A 7-*F-~ 

13 

O. 129597 

-0 . 170930E-O 7 

0. 75CC4 r E- A 1 

-'*.11926? 

0.39313'' E- r 2 

- r . 1 21 8*3 E- r 1 

- r o 

291 2 1 0 

' . ~ A 2 ? -J A 


9 

l.nooopp 





10 

P. 281 886 

l.f ccrco 




11 

r. 121291 

C. 287180 

l.OC-ornr 



12 

-P. 280595 

P. 317122 

-P.11612C 

l,rpr?in 


1 i 

-r. 21671 7 

-f .351822 

0.8C780C E-08 

-0.151453 

! # "orr ^2 


*01* UNIT05. E i IF - 


REC= FIL = 


2 
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CREDUC PROGRAM 


In optimum -seeking experimentation involving many independent variables, quadra- 
tic response surfaces are often used. A development of the most important aspect of 
the design and analysis of response surface experiments can be found in Davies (ref. 9) 
and Box and Hunter (ref. 10). A discussion of the interpretation of a quadratic surface 
fitted to a large experiment is given in reference 11. 

The general form of a quadratic surface is given by 


y = b Q + b’X + X'BX + e 



The analysis of such an equation is simplified by two calculations: (1) the calcula- 
tion of the stationary point of the surface and (2) a transformation of axes to new inde- 
pendent variables which changes the B matrix to a diagonal matrix. 

The stationary point of the surface is the solution X g to 

iy=o i = i, P 

3x i 

The transformation of the axes is given by the computation of the orthogonal ma- 
trix P which reduces B to a diagonal matrix; that is, 
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P'BP = 


h °\ 

V° V 

The A. are the eigenvalues of B, The new variables are given by 

Z = P’X 

where P is the matrix whose columns are the eigenvectors of B. If two successive 
transformations of the variables are made as follows: 

W = X - X c 

b 


Z = P'W 


then equation (26) becomes 


y ~ y s " X 1 Z 1 + X 2 Z 2 + 


+ Vp 


(27) 


where 


y s - b o + b ' x s + X ’ S BX S < 28 > 

From examination of equation (27), some general conclusions can be drawn concerning 
the attainment of a maximized response. For example, consider just two of the possible 
results. 

(1) Suppose all the A. < 0 and X is near or in the region of X at which the ex- 

periments were performed. Then clearly any deviation of Z from Z = 0 will decrease 

the response. Thus Z = 0 (or equivalently X = X ) is a maximum and is the combina- 

s 

tion of independent variables the experimenter seeks. 

(2) Some A. < 0 and some A, > 0, and X is close to the region of experimenta- 

X X s 

tion. Thus X represents what is sometimes called a saddlepoint. Moving in some 
directions will cause a decrease in y and moving in other directions will cause an in- 
crease in y. Thus the experimenter could move from X in the direction that corre- 
sponds to the direction of the Z which has the largest positive coefficient in equa- 
tion (27). This will increase y most rapidly from the value of y . 

b 
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ONI CAL REDUC 


4 75 
• . 1 ° 
.4 7,7 
274 


7 3 0 
fe 19, 
18 8. 
34 


FORTRAN STATEMENT 
ROBLEpJ 


7 42 
2 8,91 

8 9 9 


r 0, 

[ 6 , 8,8 

34 


9.248688 




,2.3 6 7 
014 9| 


V 


kSA-C-Bto IKtV 9-U-W 




Figure 8. - Sample CREDUC input 


The INPUT (Sample input is shown in fig. 8) is as follows: 

(1) Identification: One card, all 80 columns. 

(2) The number of independent variations JFAC. (14) (JFAC must be less than or 

equal to 15) 

(3) Variable Format. One card, all 80 columns. 

(4) bg (according to FORMAT in item (3)). 

(5) b 1? . . . , b JFAC (according to format in item (3)). 

(6) b lfl 

b l, 2 
b l, 3 

b l, JFAC 


b 2, 2 

b 2, 3 b 3, 3 

b 2, JFAC * • ' b JFAC, JFAC 
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One READ statement for each line, according to the format in item (3). As an ex- 
ample, consider the following estimated equation: 

y = b Q + + b 2 X 2 + b g X 3 + b 4 X 4 


+ b nX 2 




+b 12 X l X 2 

+b 22 X 2 



+ b 13 X l X 3 

+b 23 X 2 X 3 

-*33 X 3 


+b 14 X l X 4 

+b 24 X 2 X 4 

+b 34 X 3 X 4 

+b 44 X 4 

+b 15 X l X 5 

+b 25 X 2 X 5 

+b 35 X 3 X 5 

+b 45 X 4 X 5 


y = 147. 2686 - 8. 989120 Xj - 6. 817975 X 2 - 14. 60964 X g + 9. 248688 X 4 + 14. 19698 X 5 


+1.805520 Xi 


-2.126719 X^2 -1. 475730 X 2 

-7. 314218 X X X 3 -16. 10219 X 2 X 3 +2. 274270 X 3 

+ 1. 310780 XjX 4 -2. 477188 X g X 4 -1. 289688 XgX 4 +2. 236770 X 4 

+2.664464 X^5 +1. 827434 XgXg +2. 389934 XgXg +6. 014934 X 4 X 5 


-12. 58198 Xt 


A canonical reduction of this is given in the sample output. The listing of the main pro- 
gram is supplied. The subprograms TRIANG, DGELG, EIGEN, and RECT are required. 
TRIANG and RECT are as in NEWRAP. DGELG and EIGEN are the double precision 
general linear equation and eigenvalue routines from reference 12. 


SI8FTC CRFDJC 


999 


01 MENS I ON BL ! 15 ) » BS ( 105 I , BL SAVE ( 1 5 ) , BSAVE ( 105 ) , FMT ( 14 ) , XI N (225) 
DOUBLE PRECISION BL, BS, BLSAVE, BSAVE, XI N. BZERO, YS 
DIMENSION SBL( 15) ,SBS (225) ,FMTTRI ( 14) 

DATA (FMTTRI C I ) ,1=1,4) /6H(5H R0.6HW I 5, 2 , 6HX, I 801 , 6H5. 6 ) ) / 

READ (5,1001) FMT 
WRITE! 5, 1002) FMT 
READ! 5 , 1003 ) JFAC 
WRITE! 6* 1004 ) JFAC 


1 

3 

5 

7 
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REAL) ( 5*1001) FMT 
RE AD ( 5 « FMT ) BZERO 
WRITE( 6,1005 ) BZERO 
REA0(5,FMT) (BLt I ) , I = 1 « J FAC ) 

I E = 0 

DO 5 1=1, J F AC 
IE=IE*I 
I S= I E- 1 +1 

RE AD l 5 , FHT ) ( BS J K ) , K= I S , IE ) 

DO 4 <= IS, IE 
BSAVE ( K ) =BS l K ) 

4 SBS ( K) = SNGLIBS(K) ) 

BLS4VE ( I I =BL ( I ) 

SBL ( I ) = SNGL (BL t I ) ) 

BL ( I ) = -BL ( l ) 

5 CONTINUE 

LENGTH = JFAC*t JFAC+1 )/2 
WRITE(6,1006) 

HR I TE( 6,1007 ) { SBL ( I ) , I = 1,JFAC> 
WRITE(6,1 008 ) 

CALL TRIANG(SBS,JFAC,8,FMTTRI) 

I J = 1 

X I N { I ) = BStn*2,0D0 
DO 50 I =2 , JF AC 
11= 1-1 
I I K = I 

IIJ = JF AC* 1 1 
DO 40 J = 1,U 
IJ= IJ*1 

BSAVE! IJ)=0.50D0*BS( I J) 

X I N ( 1 1 K ) - BS ( IJ) 

I I K= 1 1 K + J FAC 
I I J= I I J*1 
XINtll J)= BS ( IJ) 

40 CONTINUE 
1 1 J = 1 1 J + l 
IJ= IJ + 1 

XINIII J)=2.0D0*BS( IJ ) 

50 CONTINUE 

EPS = 1 • OE— 10 

CALL OGELG! BL.XIN, JFAC, 1 ,EPS, IER ) 
IF(IER.NE.O) WRITEI6, 1009) IER 
WRITEI6.1010) 

WRITE! 6, 1007 ) ( BL ( I ) , I = 1 , J F A C ) 

I J = 0 

YS= BZERO 

DO 153 1=1, JFAC 

YS= Y S * BL ( I ) * BL S AVE ( I ) 

DO 143 J=1.I 
IJ= IJ*1 

140 YS= YS * BL! I )*BL( J)*BS (IJ) 

150 CONTINUE 

WRITE! 6,1011) YS 

C 

I J = 0 

00 163 L=l, LENGTH 

1 J= I J*- 1 

160 SBS(IJ) = SNGL(BSAVE(I J) ) 

CALL EIGENIBSAVE.XIN, JFAC.O) 

I J = D 

DO 203 1=1, JFAC 
I J= I J+ I 

200 BL ( I ) = BSAVE(IJ) 

WRITE! 6,1012) 

WRI TEC 6,1007) ( BL ( I ) , I = 1 , JF AC ) 

WRI TE! 6,1013 ) 

J J = JFAC* J FAC 
DO 213 1=1, JJ 


8 

10 

11 

12 


24 


46 

47 

54 

55 


84 

86 

87 

88 


110 


120 


130 

131 
138 


89 



148 


210 SBSII ) = SNGL (XINl I ) ) 

CALL RECTI JFAC.JFAC, JFAC.JFAC, SBS.FMTTRI ) 

GO TO 999 

C *«■ * »*#**#**********#-********#***#**#**********#**##*#***##**#*#*#*##* 

c 

c 

1001 FORMAT ( 13 A6 t A2 ) 

1002 FORMAT ( 1H1 13A6.A2) 

1003 FORMAT t 1 4 ) 

1004 FORMAT ( 1H 14) 

1005 FORMAT (24HKTHE COEFFS ARE BZERO G16.8) 

1006 FORMAT (24HK LINEAR ) 

1007 FORMAT ( 1H 8G16.8) 

1009 FORMAT (24HK SECOND ORDER ) 

1009 FORMAT (51HKTHE SOLUTION FOR STATIONARY POINT MAY BE INCORRECT / 

X 6H I ER= 13) 

1010 FORMAT (26HKTHE STATIONARY POINT IS ) 

1011 FORMAT (20HKTHE VALUE OF YS IS G16.8) 

1012 FORMAT ( 14HKE IGENVALUES ) 

1013 FORMAT ( 14HKt IGENVECTORS ) 

END 


SAMPLE CANONICAL REDUCTION PROBLEM 
5 

THE COEFFS ARE BZERO 147.268600 

LINEAR 

-8.98912001 -6.81797498 -14.6096400 9.24868798 16.1969800 



SECOND ORDER 




ROW 

1 1.805520 




ROW 

2 -2.126719 

-1.475730 



ROW 

3 -7 . 3 i 4218 

-16. 10219 

2.274270 


ROW 

4 1.3x0780 

-2.477188 

-1.289688 

2.236770 

ROW 

5 2.664464 

1.827434 

2.389934 

6.014934 


THE STATIONARY POINT IS 

3.05811524 -1.96462563 0.12068067 -3.89057544 -0. 93 71 1 06 8 E-0 1 

THE VALUE OF YS IS 120.589275 
EIGENVALUES 

9.35986030 3.75144443 1.20220299 -8.09779024 -13.9568675 


EIGENVECTORS 


ROW 

1 

-0.3001C8 

0.519711 

0. 750815 

0.238709 

-0.138316 

ROW 

2 

-0 • 549544 

-0.275150 

-0.293631 

0.690293 

-0.244071 

ROW 

3 

0.779700 

0.809764E-02 

0. 804860E-01 

0.577909 

-0. 227038 

ROW 

4 

—0.21 9489E-02 

0.790049 

-0.582318 

-0. 22882 7E -02 

-0.191618 

ROW 

5 

0. 1C5615E-02 

0.173062 

-C. 669725E-01 

0.364045 

0.912707 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, January 25, 1972, 

132-80. 
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APPENDIX - BORROWED RO UTINES 

Some of the routines used in the programs were taken from the literature. Both 
INVXTX and TRIANG are by Webb and Galley (ref. 13), and EIGEN and HIST are from 
the IBM programmer's manual (ref. 12), as are DGELG and the double precision version 
of EIGEN. 

Listing of INVXTX and TRIANG are given here, as follows: 


SIBFTC TR I AN X 

SUBROJT I NE TRIANG(A,B,NN,NKOL, FORMAT, II) I 

DIMENSION FORMAT ( 1 ) 2 

DIMENSION A(l) ,B( I) 3 

DOUBLE PRECISION B 4 

1 FORMAT ( 1 HI ) 5 

3 FORMAT ( 1H /1H /1H ) 6 

COMMON/SMALL/DUMI15) .LIST 7 

N = NN 8 

NCOL = NKOL 9 

KLUMPS = N/NCOL 10 

C 11 

KEEPTR =0 12 

K 1 = 1 13 

K2 = NCOL - 1 1* 

K3 = NCOL 15 

IF (KLUMPS .EQ. 0) GO TO 120 16 

C II 

DO 90 KLUMP=1 .KLUMPS 18 

ITRi = KEEPTR l 9 

I = -l 20 

ILO = (KLUMP-IMNCOL + ITRI + 1 21 

DO 30 K=K 1 , K2 22 

1=1+1 23 

ITRI = ITRI + K - 1 24 

ILO □ ILO + K - 1 25 

I H I = ILO + I 26 

GO TO (26,28) , I I 27 

26 WRITE! LIST, FORMAT ) K , ( A ( J ) , J = I L 0, I H I ) 28 

GO TO 30 29 

28 WRITE! LIST, FORMAT) K , ( B ( J ) , J = I LO, I H I ) 30 

30 CONTINJE 31 

KEEPTR = ITRI + K2 32 

DO 60 K=K3 , N 33 

ITRI = ITRI + K - 1 34 

ILO = ILO + K - 1 35 

I H I = ILO + NCOL - 1 36 

GO T0(56,58),II 37 

56 WRI TE( LIST, FORMAT) K , ( A ( J ) , J = I LO , I H I ) 38 

GO TO 60 39 

58 WRI TEtLIST, FORMAT) K , ( B ( J ) , J = I L 0, I HI ) 40 

60 CONTINUE 41 

K1 = <1 + NCOL 42 

K2 = <2 + NCOL 43 

K3 = <3 + NCOL 44 

90 WRITE(LIST,3) 45 

C 46 
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fsoo oo o o nooonoo 


120 ITRl □ KEEP TR 47 

IF (Kl .GT. N) GO TO 180 48 

I = -1 49 

ILO = KLUMPS*NC0L + ITRl *1 50 

DO 153 K = K 1 » N 51 

1=1+1 52 

ITRl = ITRl + K - 1 53 

ILO = ILO + K - 1 54 

I HI = ILO + I 55 

GO TO (146,1481.11 56 

146 WRITE(LIST, FORMAT) K , ( A ( J I . J = I LO, I H IJ 57 

GO TO 150 58 

148 WRITEI LIST, FORMAT) K * ( B( J ) * J = ILO* I HI ) 59 

150 CONTINUE 60 

C 61 

180 RETURN! 62 

END 63 


SIBFTC INVXXX 

SUBROJ TINE I NVXTX ( A, NN, D, FACT) 1 

2 

ASSUMES THE MATRIX A IS SYMMETRIC AND POSITIVE DEFINITE, AND ONLY 3 

THE JPPER TRIANGLE IS STORED AS A ONE-DIMENSIONAL ARRAY IN THE 4 

ORDER A(l.l), A ( 1 , 2 ) , A(2,2), Ad, 3), A(2,3>, A<3,3), .... A ( N , N ) . 5 

NN IS THE ORDER N OF THE INPUT MATRIX A. 6 

D IS (ON EXIT) THE DETERMINANT OF A, DIVIDED BY FACTOR**NN. 7 

8 

DIMENSION A ( 1 ) 9 

DOUBLE PRECISION A.PV.F 10 

N = NN 11 

ITRl =0 12 

DO 145 K=1,N 13 

14 

ITRl = ITR1+K-1 15 

KP1 = K + l 16 

KM1 = K— 1 17 

KK = ITR1+K 18 

PV = 1 . ODO/ A ( KK ) 19 

20 

I TR2 =0 21 

IF (K-l) 150,80,50 22 

23 

REDUCE TOP PART OF TRIANGLE, LEFT OF PIVOTAL COLUMN 24 

50 DO 60 J = 1 » KM 1 25 

I TR2 = ITR2+J-1 26 

KJ = ITR1+J 27 

F = A ( K J ) *PV 28 

DO 60 I=1,J 29 

IJ = ITR2+I 30 

IK = ITRl + I 31 

60 A ( I J ) = A ( I J ) + A ( I K ) * F 32 

33 

IF (K-N) 70,120,150 34 

35 

REDUCE REST OF TRIANGLE, RIGHT OF PIVOTAL COLUMN 36 

70 I TR2 = ITRi 37 
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80 DO 113 J=KP1,N 38 

I TR3 = ITR1 39 

I TR2 = ITR2+J-1 40 

KJ = ITR2+K *1 

F = AtKJMPV 42 

DO 103 1=1, J 43 

IF ( I - K ) 90,100,95 44 

90 IJ = ITR2+I 45 

IK = I TR1 + I 46 

At I J) = AC IJ) - At IK )*F 47 

GO TO 100 48 

95 IJ = I TR2 + I 49 

ITR3 = ITR3 +1-1 50 

IK = I TR3 + K 51 

At I J ) = At I J ) - At IK )*F 52 

100 CONTINUE 53 

110 CONTINUE 54 

55 

DIVIDE PIVOTAL ROW-COLUMN BY PIVOT, INCLUDING APPROPRIATE SIGNS 56 

120 I TR 2 = I TR1 57 

DO 143 1=1, N 58 

IF (I-K) 125,130,135 59 

125 IK = ITR1+I 60 

A ( I K ) = -At IK) *PV 61 

GO TO 140 62 

C tREPLACE PIVOT BY RECIPROCAL) 63 

130 A ( KK ) = PV 64 

GO TO 140 65 

135 I TR2 = ITR2+I-1 66 

K I = ITR2+K 67 

A t K I ) = A t K I ) *PV 68 

140 CONTINUE 69 

C 70 

145 CONTINUE 71 

C 72 

150 RETURN 73 

END 74 
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